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Theory  and  Simulation 

A.  Oblique  electron  Bernstein  wave  investigations.  Simulations  of  a  plasma 
driven  by  these  waves  are  continued.  The  thermodynamics  of  the  simulations 
are  checked  using  the  fluctuation-dissipation  theorem.  The  behavior  of  the 
waves  is  compared  to  linear  dispersion  theory,  and  the  onsets  of  nonlinear 
effects  are  des''»'ibed. 

B.  oiMiulation  of  the  effect  of  large  amplitude  RF  waves  on  the  interchange 
instability.  The  theoretical  groundwork  has  been  laid  for  both  the  theory 
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and  simulation  of  the  stability  of  the  interchange  mode  in  the  presence  of  an 
RF  wave.  Physically  intuitive  properties  of  the  simulation  ^uations  have 
been  found  and  estimates  of  RF  wave  quantities  have  been  derived. 

C.  One-beam  Alfven  ion-cyclotron  instabilities  of  multibeam  ion  distribution 
functions"  It  is  shown  that  for  large  wave  numbers  k,  the  Alfven  ion-cyclo¬ 
tron  instability  in  the  presence  of  a  multSpl e-beam  distribution  function  is 
dominated  by  the  response  of  a  single  resonant  beam. 

D.  Linear  mode  coupling  in  simulations  of  the  Alfven  ion-cyclotron  instability. 
Simulations  of  the  Alfvdn  ion-cyclotron  instability  show  that  modes  with 

k  >  strongly  coupled  to  each  other.  This  coupling  arises  from  | 

spatial  non-uniformity  due  to  the  particle  discreteness.  A  linear  theory  of 
the  coupling  is  derived. 

E.  A  model  of  the  plasma-sheath  region  including  ion  reflection.  Present  simu- 
lation  studies  indicate  that  increasing  ion  reflection  increases  the  electroH 
static  potential  drop  across  the  sheath  region.  A  theory  is  derived  for  the 
dependence  of  the  potential  on  the  reflection  coefficient  for  cold  and  warm 
ions.  Regimes  where  this  theory  is  not  applicable  are  described  qualita¬ 
tively. 

F.  Planar  magnetron  discharges.  We  are  using  two  approaches  to  understand 
crossed-neld  (magnetron)  discharges:  experiments  and  computer  simulations. 
We  are  proceeding  with  the  design  and  construction  of  a  planar  discharge  of 
potentials,  fields,  densities  and  velocities  as  a  function  of  a  variety  of 
parameters.  We  are  also  using  plasma  computer  simulation  for  a  similar  con¬ 
figuration  in  order  to  measure  the  same  qualities. 

G.  Particle  simulations  of  the  low-alpha  Pierce  diode.  Particle  simulations 
are  used  to  study  the  evolution  of  small  initial  perturbations  about  the 
classical  Pierce  diode  equilibrium  in  the  parameter  range  0  <  a  =  mpL/vQ  <  3tt. 
Linear  properties  recovered  show  quantitative  agreement  with  theoretical 
predictions. 

H.  Theory  and  simulationof  ion-acoustic  double  layers.  The  formation  of  double 


layers  is  investigated  under  the  conditions  m./m  = 

I  C 


100  and  T 


is  found  that  weak  ion-acoustic  double  layers  can  form  even  in  short  systems 
(80  Xg)  and  with  relatively  small  electron  drift  (v  =  0.45  v^j^  =  4.5  c  ). 
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ELECTRONICS  RESEARCH  LABORATORY 

University  of  California 
Berkeley,  California  94720 


SECTION  I;  PLASMA  THEORY  AND  SIMULATION 


A.  Oblique  Electron  Bernstein  Wave  Investigations 

Wm .  S  .  Lawson 


Much  progress  has  been  made  in  the  investigations  of  electron  Bernstein 
waves  begun  last  year  using  a  modiflcation  of  the  code  PDWl.  Some  results  of 
interest  are  here  detailed. 

BWADR  —  linear  dispersion  relation  code 

The  BWADR  code  (due  to  Jens-Peter  Lynov),  was  written  to  compute  the 
complex  dispersion  curves  for  oblique  Bernstein  waves  The  resulting  curve 
show  some  interesting  features,  and  are  worthy  of  discussion.  The  parameters 
chosen  are  Op/Ue=2  and  i>=80°.  The  real  and  imaginary  parts  of  A:  are  plotted 
for  real  o. 

Figure  la  shows  the  frequency  region  between  0  and  Uf  There  is  no  mode 
in  this  region  for  normal  Bernstein  waves  (i5=90°).  The  mode  is,  however,  deriv¬ 
able  from  cold  MHD  theory.  It  has  the  dispersion  relation  (which  includes  the 
upper  hybrid  mode) 

u*  -  (up-hUe)u^  +  w^t>|cos®i>  =  0 

These  modes  are  plotted  versus  tJ  in  figure  lb.  The  lower  mode  degenerates  to  a 
static  mode  at  15=90”,  and  to  an  E=0  mode  at  i5=0“.  This  mode  may  well  be 
unobservable  for  normal  plasma  parameters  when  electromagnetic  effects  are 
taken  into  account. 

Figure  2  shows  the  frequency  region  between  Ug  and  2Uj.  The  normal 
Bernstein  dispersion  curve  is  shown  for  comparison.  For  a  time,  the  two  are 
very  close,  but  for  larger  k  values  the  damping  rate  becomes  large,  and  the  real 
part  of  k  shifts  drastically.  This  branch  of  the  wave  was  used  for  most  of  the 
simulations. 

Figure  3a  shows  the  next  frequency  interval,  that  between  Zug  and  ZUg. 
This  interval  shows  fuio  modes  of  importance.  The  structure  of  the  modes  is 
especially  interesting  in  comparison  to  the  normal  Bernstein  wave  curve  (figure 
3b).  The  forward  and  backward  portions  of  the  curve  for  normal  waves  splits 
into  two  branches,  each  of  which  is  purely  forward  or  purely  backward.  Both 
modes  have  large  damping  rates  above  the  maximum  in  u  for  normal  waves,  and 
the  backward  wave  has  large  damping  below  that  maximum.  The  first  part  of 
the  forward  wave,  however,  has  very  little  damping.  This  wave  corresponds  to 
the  upper  hybrid  wave  of  MHD,  and  is  the  dominant  wave  in  the  noise  spectrum. 

In  figure  4.  it  can  be  seen  that  the  higher  harmonics  are  rapidly  becoming 
unimportant  because  of  very  strong  damping.  Higher  harmonics  are  so  strongly 
damped,  that  they  can  be  for  all  intents  and  purposes  ignored. 

Fluctuation  Dissipation  Cheek 

One  problem  which  became  apparent  early,  was  the  high  noise  level  in  the 
simulations,  even  for  large  numbers  of  particles.  This  suggested  a  method  of 
checking  the  operation  of  the  code,  namely  to  compare  the  noise  level  with  the 
level  predicted  by  the  fluctuation-  dissipation  theorem.  This  theorem  states 
that 
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A  similar  method  was  used  successfully  by  Kamimura,  Wagner,  and  Dawson  [l]  in 
two-dimensional  simulations  of  perpendicular  Bernstein  waves. 

The  post-processor  ZED  (with  some  modification)  was  used  to  compute  the 
fluctuation  spectrum  of  the  electric  field  at  a  point  one  quarter  of  the  way 
across  the  system  from  the  source.  The  result  is  seen  in  figure  5.  Plotted  are 
the  theoretical  and  observed  values  of  <\  E{d)\^> /  ZkT  versus  frequency  (note 
that  for  a  given  Up  andu,.  T  is  a  function  of  the  particle  number  density). 

The  numerical  integration  of  the  theoretical  spectrum  deserves  some  comment, 
as  it  is  not  simple.  The  integrand  has  resonances  in  it,  some  of  them  very  nar¬ 
row  and  strong.  These  resonances  must  be  subtracted  to  two  orders  in  the 
width  parameter  in  order  to  get  something  smooth  enough  to  integrate  numeri¬ 
cally.  Even  this  function  is  less  than  ideal,  but  the  algebra  of  subtracting 
another  order  is  much  more  difficult. 

The  resolution  is.  of  course,  much  higher  in  the  theoretical  plot,  but  the  degree 
of  agreement  can  be  seen  by  comparing  a  broad  peak  such  as  the  one  at 
«=3.2cje.  One  apparent  point  of  disagreement  is  the  appearance  of  peaks  at  Uc 
and  2ue  in  the  simulation.  The  origin  of  these  peaks,  if  they  are  indeed  real,  is 
unknown,  but  the  peaks  may  be  due  to  the  finite  size  particle  approximation. 
This  agreement  between  the  theoretical  and  observed  spectra  gives  some 
confidence  in  the  accuracy  of  the  simulation. 

To  see  the  magnitude  of  the  noise  relative  to  the  signal  in  a  typical  driven  run, 
figure  6  shows  a  spectrum  plot  of  such  a  run  at  an  excitation  amplitude  which 
just  begins  to  show  some  non-linear  effects.  All  these  simulations  used  130,000 
particles  —  a  large  number  for  a  one  dimensional  simulation.  A  heuristic  reason 
for  the  high  noise  level  is  the  large  number  of  phase  space  dimensions  (four) 
which  must  be  filled.  Just  as  increasing  the  number  of  spatial  dimensions 
increases  noise,  so  does  increasing  the  number  of  velocity  dimensions. 

linear  and  Non-linear  Amplitudes 

To  find  the  limits  of  linear  behavior,  and  the  character  of  non-linear 
behavior,  a  set  of  runs  was  made  at  one  set  of  parameters  for  varying  excitation 
amplitudes.  The  observed  spatial  damping  was  then  compared  to  linear  theory 
and  examined  for  interesting  phenomena.  These  parameters  are  11=80°. 
t)p/Ue=Z,  =  (uo  is  the  driving  frequency),  and  the  length  of  the  sys¬ 

tem  was  sixty-four  thermal  gyroradii.  The  results  were  rather  interesting,  and 
will  be  discussed,  but  first  the  method  of  obtaining  the  diagnostic  plots  must  be 
discussed,  as  the  diagnostic  is  not  straight-forward. 

Because  of  the  high  noise  level,  even  very  non-linear  signals  would  be  unobserv¬ 
able  unless  the  noise  were  filtered  out  somehow.  The  method  used  to  do  this 
was  to  wait  until  the  system  had  come  to  equilibrium,  then  take  a  running 
Fourier  transform  af  ths  driving  frequency  only.  This  results  in  plots  of  the 
real  and  imaginary  parts  of  £(x,&;o).  These  profiles  were  then  squared  and 
added,  and  plotted  on  a  log  scale,  yielding  a  plot  of  ln|£'(z,uo)r  versus  z. 
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Linear  theory  predicts  that  this  plot  should  be  a  straight  line. 

The  lowest  amplitude  excited  is  plotted  in  figure  7.  The  deviations  from  a 
■traight  line  are  small,  and  show  no  structure.  The  over-all  slope  of  the  curve  is 
a  very  good  match  to  the  straight  line  of  linear  theory.  Thus  it  is  fair  to  con¬ 
clude  that  the  system  is  behaving  in  a  linear  fashion,  and  that  linear  theory  is 
accurate  in  some  neighborhood  of  the  chosen  parameters. 

The  next  plot  (figure  8)  shows  an  excitation  amplitude  twice  that  of  the  first. 
The  decay  rate  is  still  a  good  match  with  the  linear  decay  rate,  but  beats  and  a 
plateau  have  appeared.  Since  these  beats  disappear  when  the  left  wall  is  made 
absorbing  (and  emitting)  rather  than  reflecting,  the  beats  can  be  ascribed  to  a 
ballistic  mode  which  has  traveled  through  the  exciter,  i.  e.  the  electrons  have 
been  bunched  in  the  process  of  damping  the  wave  (which  has  a  phase  velocity 
moving  toward  the  exciter). 

The  third  plot  (figure  9)  is  again  at  twice  the  amplitude  of  the  previous  one. 
Now  a  region  of  reduced  damping  near  the  exciter  becomes  evident.  This  is 
almost  certainly  evidence  of  the  onset  of  trapping.  The  wave  can  be  damped 
only  to  a  point  at  the  linear  rate,  then  a  slower  process  must  take  over. 

This  becomes  clear  in  the  next  plot  (figure  10),  which  is  at  five  times  the  previ¬ 
ous  excitation  level.  Now  the  electrons  are  incapable  of  damping  any  but  a 
small  fraction  of  the  wave  energy  (except  at  the  wall,  where  the  damping 
processes  are  presumably  of  a  different  nature).  The  wavelength  of  this  wave 
has  also  roughly  doubled,  though  this  may  be  due  to  the  heating  of  the  distribu¬ 
tion.  The  reason  for  the  heating  is  also  unclear.  It  may  be  due  to  the  passage 
of  the  wave,  but  it  seems  more  likely  that  it  is  due  to  the  exciter  itself. 

The  final  figure  (figure  11)  shows  an  excitation  amplitude  of  ten  times  the  previ¬ 
ous  one.  Now  the  damping  returns  at  almost  the  linear  rate,  but  the  process 
must  be  different.  The  most  likely  explanation  is  a  surfing  effect  which  occurs 
when  the  electric  field  is  large  compared  to  vxB  [2].  This  effect  can  result  in 
large  damping  even  for  Bernstein  waves  at  90°  to  the  magnetic  field. 

[1]  Kamimura,  Wagner,  and  Dawson.  Phys.  Fluids  81.  1167  (1978) 

[2]  T.  Katsouleas.  Ph.D.  Thesis,  1984.  UCLA 
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SIMULATION  OF  THE  EFFECT  OF  LARGE  AMPLITUDE 
RF  WAVES  ON  THE  INTERCHANGE  INSTABILITY 
—SUPPORTING  THEORY 

Niels  F.  Otani  (Dr.  J.  A.  Byers  and  Prof.  C.  K.  Birdsall) 


I.  Introduction 

Recent  experiments  at  the  Phaedrus  tandem  mirror  experiment  at  Wisconsin*  have 
provided  evidence  that  imposed  RF  waves  with  frequencies  in  the  vicinity  of  the  ion- 
cyclotron  frequency  may  have  a  strong  stabilizing  effect  on  interchange  modes  in  an  ax- 
isymmetric  plasma.  In  the  experiments,  the  stabilizing  properties  of  the  RF  wave  were 
observed  to  change  sharply  at  a  critical  frequency  near  the  ion-cyclotron  frequency,  with 
the  strong  stabilizing  effect  evident  only  above  the  critical  frequency.  It  was  suggested* 
that  the  stabilizing  effect  may  be  due  a  ponderomotive  force  produced  by  the  RF  wave^ 
which  was  counteracting  the  effective  gravity  driving  the  interch2uige  instability. 

In  view  of  recent  results,  we  have  attempted  to  simulate  some  of  the  effects  which  may 
be  contributing  to  the  stabilization  of  the  interchange  mode  in  the  Wisconsin  experiment. 
The  simulation  code  used  is  a  modification  of  the  2  1/2-d  quasineutral  Darwin  code  PEPSI 
described  in  a  previous  report.^ 

The  code  simulates  the  x-y  plane  perpendicular  to  the  magnetic  field  B(x,  y)i.  The 
RF  wave  propagation  vector  ko(z),  the  equilibrium  density  Vno(x)  and  an  external  grav¬ 
itational  field  g(z]  are  all  initialized  to  point  in  the  x-direction,  while  a  small  interchange 
perturbation  is  initialized  with  a  sinusoidal  dependence  in  the  y-direction.  Full  ion  dy¬ 
namics  are  included  in  the  code,  while  rij  =  n*  and  u*  =  cE  x  B/B*  are  assumed  for  the 
electron  density  and  fluid  velocity  respectively.  A  complete  description  of  the  code  along 
with  simulation  results  will  appear  later;  however,  individual  results  will  be  quoted  here  as 
needed. 

In  this  report,  we  present  supporting  theory  for  the  simulation.  Appropriate  fluid 
equations  for  our  simulations  will  be  derived  and  some  of  their  properties  demonstrated. 
Characteristics  of  the  RF  wave  to  second  order  in  the  wave  amplitude  will  be  described. 
Equations  describing  the  linear  behavior  of  the  interchange  mode  in  the  presence  of  the  RF 
wave  are  stated. 


II.  The  Equations 

The  algorithm  used  by  our  simulations  may  be  represented  by  the  fluid  equations  : 


du  _  e  u  X  B  \ 

—  -t-uVu  =  — (E-t- - )  g, 

at  m\  c  / 


dn  T7  r  \ 

=  -V  (-»). 


(Id) 


/V  xB 
\  4iren 


X  B, 


where  u  is  the  ton  fluid  velocity,  n  =  =  n*  is  the  density,  and  m  is  the  ion  mass.  The 

form  of  the  ion  momentum  equation  Eq.  (la)  is  appropriate  because  all  simulations  to  date 
have  been  initialized  with  cold  ions.  Equation  (Id)  comes  from  the  zero- pressure,  inertialess 
electron  momentum  equation 


0  =  E  + 
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(2) 


where  for  u*.  Ampere’s  Law  has  been  substituted  neglecting  the  displacement  current. 

Since  the  simulation  takes  place  in  the  x-y  plane,  the  magnetic  field  need  only  have 
a  2-component:  B  =  B(x,y,t)i.  Noting  also  that  E  only  appears  in  the  combination 
E  u  X  B/c,  Eqs.  (1)  may  be  simplified  to 
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where  u  =  (ux,Uy,0).  Equations  (3)  form  the  basis  for.  all  subsequent  analysis  in  this 
report. 


III.  Properties 

Several  properties  may  be  derived  for  the  system  described  by  Eqs.  (3).  A  number  of 
these  properties  are  particularly  appealing  because  they  have  a  simple  physical  interpreta¬ 
tion. 

From  Ampere’s  Law,  the  electron  fluid  velocity  may  be  expressed  as 

u«  =  u  -  -^—VB  X  i,  (4) 

4jren 

from  which  immediately  follows  J  •  VB  =  0  and  V  •  J  =  0,  as  is  also  true  of  MHD  theory. 
Thus  the  total  current  flows  along  constant-B  contours  and  is  divergence-free. 
Substituting  Eq.  (4)  into  Eq.  (3b)  yields 

=  -V  •  (nux),  (5) 

the  electron  continuity  equation,  while  substitution  into  Eq.  (3c)  produces 


demonstrating  conservation  of  magnetic  flux  /  Bdxdy.  Equations  (3a),  (4),  (5),  and  (6) 
are  an  alternative  and  equivalent  set  of  equations  to  Eqs.  (3).  The  similarity  of  Eqs.  (5) 
and  (6)  may  be  exploited  to  yield 


(7) 


i.e.,  the  ratio  Bfn  is  constant  in  the  electron  fluid  frame.  Thus,  the  magnetic  field  may  be 
thought  of  as  being  frozen  into  the  electrons. 

An  energy  conservation  theorem  also  exists  for  this  system.  Taking  the  dot  product  of 
both  sides  of  Eq.  (3a)  with  mnu  produces 


1  dn}  1  , 

-mn——  +  -mnu  -  V« 

2  at  2 


-  ~uB  ■  VB  +  mnu  •  g. 
4x 


(8) 


Ampere’s  Law  (Eq.  (4))  may  be  used  to  show 


u5  •  VB  =  u*fl  •  VB.  (9) 

Multiplying  both  sides  of  Eq.  (6)  by  B  yields 

u.BVB=  ~(B=‘)  +  V(ueB^).  (10) 

Using  ion  continuity  on  the  left  side  of  Eq.  (8)  and  substituting  Eqs.  (9)  and  (10),  we  obtain 

—  ( -nmti^  -I-  nm<t>g  I  +  V  •  ( -nmw'u  +  — u*  +  nm(^,u  1=0,  (11) 

where  g  =  —V0g.  Thus  the  local  total  energy  density  define<.  as  nmu*/2  +  B^fSir  +  nmd>g 
will  be  conserved  in  any  region  on  the  x-y  plane  up  to  a  flux  on  the  region  boundary. 


IV.  The  RF  Wave 


Estimates  of  the  magnitudes  of  various  quantities  associated  with  the  RF  wave  will  be 
required  for  the  stability  analysis  of  the  interchange  mode  when  the  RF  wave  is  present. 
Since  effects  second  order  in  the  RF  amplitude  are  expected  to  contribute  the  main  stabi¬ 
lizing  influence,  calculations  to  that  order  were  performed.  We  use  the  equations 


dux  dux 

+  «x- 


dt 


dx 


B  dB 
4irmn  dx 


+  9. 


dUy  dUy 
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since  the  RF  wave  has  been  assumed  to  depend  only  on  z. 


(12a) 

(12b) 

(12c) 
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Zero  order.  The  equilibrium  is  assumed  to  be  time-independent  [djdt  =  0)  with 
Ui  =  Uxo  =  0-  Equation  (12a)  then  implies  a  depression  must  exist  in  the  magnetic  field 
structure  due  to  the  gravity-induced  current: 


±(?l\  = 

dx\%it) 


while  Eq.  (12b)  allows  to  be  arbitrary.  If  we  also  require  the  electric  field  Eq  to 

be  zero  in  equilibrium,  Eq.  (la)  forces  u,o  =  —mcqjeBo,  the  gravitational  drift  velocity. 
Equations  (12c)  and  (12d)  yield  no  further  information. 

First  order.  Equations  (12a)  through  (12c)  may  be  linearized  about  the  equilibrium 
to  yield  a  modified  wave  equation: 

where  u^i  is  the  first  order  component  of  u*,  =  Bo/(4nmno),  Lg^  =  (1/ Bo)(dBoldx), 

and  =  (l/no)(3no/3z).  Equation  (12d)  may  be  linearized  to  yield  an  equation  for  Uyi 
which  evidently  plays  no  role  in  the  wave  mechanics. 

For  the  remainder  of  this  section,  a  uniform  equilibrium  (ff  =  0,  dnojdx  —  0)  will  be 
assumed,  because  as  stated  earlier,  only  estimates  of  the  wave  quantities  will  be  required. 

In  a  uniform  system  Lg*  =  =0;  thus  Eq.  (14)  becomes  simply  the  wave  equation 

for  the  propagation  of  compressional  Alfven  waves.  If  we  assume  a  standing  wave  of  a 
given  frequency  uq  (this  is  the  type  of  wave  initially  loaded  in  the  simulation),  we  find  the 
following  relative  amplitudes  and  spatial  and  temporal  phases  among  the  wave  components: 

Si  =Bcoskoxcosu;ot,  (lha) 


Til  =no-^  coskoxcoauot,  (ISb) 
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where  fco  =  ^oI^a- 

Second  order.  When  Eq.s  (12a)  and  (12d)  are  expanded  to  second  order  and  the 
uniform  equilibrium  assumption  is  again  made,  the  result  is 


dt^  I  So 


where  expressions  from  Eqs.  (15)  have  been  used.  This  has  a  particular  secular  solution 


=  i/'AV 
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which  has  been  observed  in  simulation.  The  d.c.  portion  of  the  solution  is 


-Hi.)' 


cos2A:oz. 


The  equation  corresponding  to  Eq.  (16)  may  be  derived  for  u^a;  it  has  a  secular  solution 
whose  d.c.  part  is 

(urf}=0.  (19) 

The  expression  for  (Ba)  in  Eq.  (18)  provides  some  insight  into  the  nature  of  the  second 
order  current.  Since  (u^o)  =  i^eyo)  =  («yi)  =  0,  we  have 

4ir  4jreno.  .  iire  ,  . 

- r -  =  - {jyil  —  - {“v2  ~  “eya) - (^lUeyl)-  (20) 

OX  c  c  c 

The  last  term  in  Eq.  (20)  may  be  evaluated  using  Eqs.  (15): 


1  /  B  \  ^ 

(niUeyi)  = -koBol—j  sin2*:ox, 


which  we  observe  from  Eq.  (18)  is  equal  to  -d{Bi)/dx.  Thus 

(ttya  —  Ueya)  =  0,  (22) 

which  is  somewhat  of  a  surprise.  We  might  have  expected  the  second  order  current  would 
be  provided  by  a  ponderomotive  ion  drift  of  the  form 

(«,,)  =  (23) 

where  Fp  is  the  ponderomotive  force  due  to  the  RP  wave;  however  if  this  is  the  case,  it  is 
canceled  by  the  second  order  electron  drift,  and  instead  (riiUeyi)  must  be  considered  as  the 
important  second-order  current  term.  It  is  not  clear  whether  this  second-order  mechanism 
will  be  of  the  type  necessary  to  reverse  the  destabilizing  gravitational  current. 


V.  Stability  of  the  Interchange  Wave 

By  considering  the  system  containing  the  RF  wave  to  be  a  modified  equilibrium,  it  is 
possible  to  investigate  the  stability  properties  of  the  interchange  mode  in  the  presence  of  the 
wave.  Assuming  the  interchange  perturbation  to  be  of  the  form  dn{x,y,t)  =  6n{x,t)e'^^ 
(similarly  for  6u  and  6B),  we  obtain  to  first  order 
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All  equilibrium  quantities  in  Eqs.  (24)  contain  the  RF  wave  and  therefore  have  both  space 
and  time  dependence.  Here  (*;„  =  eB/mc  and  G  ~  g  —  dxizjdt  —  Uzduxidx.  Equation  (24) 
should  be  solvable  in  the  local  approximation  using  a  method  analogous  to  that  used  to 
solve  Mathieu's  Equation  for  multiply  periodic  solutions;  this  calculation  is  still  in  progress. 

In  the  case  of  a  vanishingly  small  RF  wave,  the  assumptions  k'^v\  <  and  tj  <S  ujc, 
in  the  local  aproximation  yield  the  usual  interchange  dispersion  relation 


-  2  ,  ^9  -  ,  9  ^^0  ^ci  _  n 

U)  +•  - UJ  + - — - =  0 

Wc»  no  ox  (jJci 


(25) 


where  Q  =  uj  —  kuoy,  w  is  the  interchange  mode  frequency,  and  Wc,  =  ujd  +  duoyfdx  is  an 
effective  ion  cyclotron  frequency. 


VI.  Conclusions 

Theoretical  groundwork  has  been  laid  for  both  the  theory  and  simulation  of  the  stabil¬ 
ity  of  the  interchange  mode  in  the  presence  of  an  RF  wave.  Physically  intuitive  properties 
of  the  simulation  equations  have  been  found.  Estimates  of  RF  wave  quantities  have  been 
derived.  We  find  that  the  current  second-order  in  the  RF  wave  flows  due  to  an  unexpected 
mechanism  related  to  first  order  perturbations  in  the  density  and  the  electron  fluid  velocity. 
Finally  equations  linear  in  the  interchange  perturbation  have  been  derived  and  are  in  the 
process  of  being  solved. 
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C .  ONE-BEAM  ALFVfiN  ION-CYCLOTRON  INSTABILITIES 
OF  MULTIBEAM  ION  DISTRIBUTION  FUNCTIONS 


Niels  F.  Otaoi  (Prof.  C.  K.  Birdsall) 


Tests  of  the  simulation  code  TRACY  used  to  study  the  Alfven  ion-cyclotron  (AIC) 
instability  have  included  simulations  using  multibeam  ion  distributions  of  the  form, 


/(v) 
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2irv±Nk 
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^(t'X  -  -  Vx6). 

6=0 


(1) 


where  -Vj,  is  the  number  of  beams.  These  distributions  were  observed  in  simulation  to  yield 
growth  rates  very  close  to  those  expected  of  bimaxwellian  ion  distributions  having  the  same 
T,||  and  T,±  for  wavenumbers  k  on  the  principal  growth  peak.  However,  a  deviation  from 
bimaxwellian  growth  rates  was  observed  for  high  k  modes.  In  this  article,  we  demonstrate 
that  the  observed  deviation  of  the  growth  spectrum  for  high  k  may  be  explained  by  “one- 
beam”  AIC  instabilities.  That  is,  for  high  enough  k,  it  is  possible  for  the  ion  response  to 
be  dominated  by  a  single  beam  in  the  ion  distribution  function.  The  resonant  beam  has 
infinite  temperature  anisotropy  and  therefore  is  capable  of  driving  the  AIC  instability. 

From  a  form  of  the  AIC  dispersion  relation  derived  previously,*  we  easily  obtain  the 
disp*.  'ion  relation  for  the  multibeam  distribution  function  described  by  Eq.  (1); 
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We  will  be  looking  for  a  mode  associated  with  a  given  beam,  say,  the  beam  6  =  0  in 
Eq.  (1).  Accordingly,  we  rewrite  Eq.  (2)  as 
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where  wd  =  w  —  —  fl  is  the  deviation  of  the  mode  frequency  from  the  Doppler-shifted 

ion-cyclotron  frequency.  The  mode  will  be  assumed  associated  with  the  beam  6  =  0  in  the 
sense  that 

Wd  =  Wcio  +  ^  ^  as  fc  —  00,  (4) 

where  ujjq,  wdi,  and  Wdz  are  independent  of  k;  i.e.,  we  are  looking  for  roots  near  the  resonant 
frequency  of  beam  6  as  fc  — ►  oo. 

With  this  form  of  u/d.  the  dispersion  relation  Eq.  (3)  may  be  expanded  in  powers  of 
fc~‘  as 


■T-»T 


i 


I 

k 

i 


I 


k 


r '. 

r. 

r.* 


to 


P 


f>>\ 

m 


V 

V 


k 

■l 


+  0(k-^). 


(5) 


Solving  Eq.  (5)  power  by  power  in  fc~‘,  we  find  that  the  k^,  k^,  and  k°  terms  yield  respec¬ 
tively, 
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Substituting  Eq.  (6)  into  Eqs.  (7)  and  (8),  and  then  substituting  Eq.  (7)  into  Eq.  (8),  we 
find 
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Equations  (6),  (9),  and  (10)  are  thus  the  terms  in  Eq.  (4)  which  described  the  behavior  of 
individual  beam-related  modes  as  1:  — '  oo. 

For  large  enough  k,  the  u;ao  term  dominates  in  Eq.  (4)  and  we  recover  as  the  imaginary 
part  of  the  frequency  just  the  infinite-anisotropy  AIC  growth  rate  for  that  beam.  The 
expression  for  wao  is  pure  imaginary,  implying  that  the  real  part  of  the  mode  frequency  is 
just  the  Doppler-shifted  ion  cyclotron  frequency,  Q  -1-  fcv,o-  These  being  the  charzurteristics 
expected  of  the  AIC  instability  for  an  ion  distribution  consisting  of  the  one  beam,  we 
conclude  that  for  large  enough  k,  Doppler-shifted  AIC  instabilities  do  indeed  exist  for  the 
individual  beams  of  the  ion  distribution. 
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1).  LINEAR  MODE  COUPLING  IN  SIMULATIONS  OF 
THE  ALFV6N  ion-cyclotron  INSTABILITY 

Niels  F.  Otani  (Prof.  C.  K.  Birdsall) 


The  identical  linear  growth  rates  of  the  short-wavelength  (jkvyt/fl  >  1)  spatial  Fourier 
modes  observed  in  Alfven  ion-cyclotron  (AIC)  simulations  (other  than  those  involving 
multibeam  ion  distribution  functions)  suggest  that  the  various  wavelengths  are  strongly 
coupled.  Such  coupling  in  the  linear  regime  only  occurs  when  spatial  non-uniformities  are 
present  in  the  equilibrium.  In  this  case,  the  non-uniformities  arise  from  the  spatial  ion- 
particle  distribution.  These  non-uniformities  are  inevitable  in  warm  particle  simulations, 
since  even  if  particles  are  initially  lo2ided  with  regular  spacing,  the  regularity  is  immediately 
destroyed  as  the  simulation  pushes  the  particles  with  their  various  zero-order  velocities. 
There  is  one  exception;  if  each  velocity  to  be  loaded  is  represented  by  a  set  of  particles 
uniformly  distributed  throughout  the  simulation  (i.e.,  a  velocity  beam),  then  clearly  the 
uniformity  of  the  zero-order  particle  distribution  will  be  maintained.  This  explains  why  our 
multibeam  AIC  simulations  do  not  suffer  from  coupled  linear  growth  of  the  high-fc  modes. 

The  non-uniformity  of  the  spatial  particle  distribution  enters  the  equations  through  lin¬ 
earized  source  terms.  In  the  case  of  our  AIC  simulations,  the  relevant  term  is  the  linearized 
transverse  ion  current.  The  essence  of  the  linear  coupling  is  thus  easily  understood.  The 
linear  current  response  of  ions  to  a  transverse  wave  at  a  given  wavelength  will  not  be  en¬ 
tirely  at  that  wavelength  because  the  constituent  linear  particle  currents  are  not  uniformly 
distributed  in  space. 

The  details  of  this  linear  coupling  effect  may  thus  be  formulated  by  examining  the 
effect  of  the  linear  current  response  on  the  field  equations.  The  linearized  expression  for 
the  current  J  at  grid-point  J  used  in  both  simulation  codes  TRACY  and  PEPSI*  is 

Jy (0  =  E  E  E 

dt 

where  no  is  the  unperturbed  density,  Ng  is  the  number  of  the  grid-points  in  this  periodic 
system,  the  index  i  refers  to  groups  of  particles  with  identical  zero-order  values  of  z,  vm,  and 
Vx,  a  indexes  the  particles  in  each  such  group,  N  is  the  total  number  of  particles,  zj’(<)  = 
2°  +  Uyf,  u°,(<)  =  sin(fl°,  +  fit)  +  "i"  ^0  the  zero-order  trajectories, 

S(2)  is  the  shape  factor  of  each  particle,  and  zj  is  the  z-position  of  grid-point  J.  Here 
6v,,(t)  and  6zi,{t)  are  the  linear  perturbed  quantities  for  particle  a  in  group  i,  Aj  is  the 
perturbed  vector  potential  at  grid  J,  Bq  =  Bo*  is  the  unperturbed  magnetic  field  which 
is  uniform,  constant  in  time  and  parallel  to  the  grid  direction,  and  H  is  the  corresponding 

*  Note:  PEPSI  does  not  use  the  vector  potential  A.  Thus  slight  errors  arising  from 
finite  differencing  are  introduced  here.  These  will  be  ignored. 


ion-cyclotron  frequency.  Again,  ions  are  represented  as  particles  while  electrons  enter  only 
through  the  electron  E  x  Bq  current.  In  this  calculation  we  ignore  the  finite  timestep  size 
At. 

Expressions  for  X!,  5v,j(t)  and  6z„(t)  in  the  ion  terms  in  Eq.  (1)  may  be  calculated 
from  the  linearized  ion  equation  of  motion  : 

+  n  X  =  ~Y. 

>  ' 

^  E  [(I  +  A*  -  Vk(v“  •  Afc)]e‘**>5(z“(t)  -  z,)(2) 

where  we  have  rewritten  Aj(f)  =  51*  A*:(f in  which  j  is  a  grid-point  index  and  k 
assumes  the  values  k  =  2iTmlL,  m  =  -lVj/2  -t- 1, . . . ,  Ng/2  and  k  =  ki. 

Considering  first  the  perpendicular  perturbed  velocity  component  =  i^Viax  + 
i^Viay)/\/2,  we  find 

^  E  (  ~  ^  ~  ikv°^^A^(t)\SkNgexp[ik{^°  -I-  v°||t)l,  (3) 

where  the  particle  shape  factor  has  been  Fourier  transformed  S(z)  =  5Z*=_oo  and 

grid  aliasing  effects  have  been  ignored  : 

«SfciV,4v. 

3 

This  is  often  a  good  approximation  for  wavenumbers  k  of  interest,  i.e.,  wavenumbers  well 
within  the  first  Brillouin  zone  -ir/Az  <  A:  <  tt/Az;  however,  evidence  suggests  that  grid 
aliasing  is  in  fact  important  for  this  type  of  simulation.  This  will  be  amplified  in  future 
reports. 

Assume  now  that  A*  (<)  is  a  sum  of  waves, 

A+(t)  =  |du;A+Me-“‘,  (4) 

with  the  values  of  (complex)  w  included  in  the  sum  to  be  specified  later. 

Equation  (3)  is  easily  solved  for  by  multiplying  by  the  integrating  factor  c*^': 

^(e‘«‘E^Vu)  =  J]  /  <fw5fclVg(w-*vf||)A;t{w)exp(jfc{z?-f-v°||t)-*wtl,  (5a) 

J  ' 

+  e~‘"‘  f  dt'  f  duSkNgi{u  -  A:t;°|)A;J’{w)exp(iA:(z°  +  u°||f')  -  (w  -  n)<']Y  (5b) 


yielding  on  integration, 

J  ' 

~Y1  j  <^^kNg  ^  (w)  explifc(2“  +  v°||t)  -  iut]^ .  (5c) 

The  expression  for  Szi^(t)  is  obtained  by  examining  the  parallel  component  of  Eq.  (2). 
Again  ignoring  grid  aliases  and  using  Eq.  (4)  we  find  the  expression  analogous  to  Eq.  (5a): 

du)  SkNgik(\°^it)  ■  A*(w))exp(ifc(2°  +  -  tu;<|.  (6) 

Since  in  the  }  basis, 

v?,(0  •  Afc(u)  =  -^^{AjJ'(uj)exp(i(n<  +d°  ))  -  A^  (w)  exp[-i(nt  (7) 

we  find  upon  substituting  into  Eq.  (6) 


X  {AjJ'(a;)c‘®'*  exp[-i((J  -  fcu°|  -  0)t\  -  A^  (w)c“**‘*  exp[-t(ui  -  fcn°|  +  n)t)}.  (8a) 

Two  integrations  with  respect  to  t  yield 


V2 


/  .  +  ,  ^  «xp(-t(u;  -  fcv°  -  n)t]  ^  expl-i(w  -  fcu°  +  n)tj 

X  (^A,  (u;)e  ••  _  A^  (u;)c  ••  _ 


-  Afc  (u;)e‘ 


-i(u;  -  fcv°||  +  n) 


= £  [c;; + c!j<  -  e/<^ 

/  .  +  ,  ^  exp[-t(u;  -  fcv°  -  n)ti  _  o  exp(-i(u; 


(  .  +  ,  ^  .e°  «p[-.(u;  -  fcv°  -  n)ti  ^  _.,o  exp(-.(u;  -  fcv°|  +  n)‘lM 

X  (■^*<">' - (ZTTijjpnjr-  -  ^*1")' - oTTltllTw— jj- 

The  (+)-polarization  portion  of  the  Szit(t)  term  in  Eq.  (1)  may  now  be  evaluated: 

x"  ■  ^  v°.j_(t)«z<,(t)  =  ^  ^c-<c-"‘6z.,(t).  (9a) 

Here  x~  =  (x'*')*  =  x  +  »y  is  the  dual  basis  vector  to  x"*"  =  x  -  iy.  Now  assuming 
the  particles  in  each  group  »  are  uniformly  loaded  in  gyrophase  space  with  at  least  three 


particles  (the  four-spokes  distributions  used  in  our  simulations  have  four  uniformly  loaded 
particles  per  group)  then 

N. 

S=\ 

kills  the  last  term  in  Eq.  (8c)  when  substituted  into  Eq.  (9a).  Thus  we  obtain 


X~  •  +C‘,||<e 

4  *  ^ 

Finally,  the  E  x  B  term  from  Eq.  (1)  may  be  evaluated  as 

dJi,  /• 

x~  •  X  Bo  =  -  ^  j  dui  ujBqA^ (u)  exp(ikzj  -  iwt).  (10) 

Substituting  Eqs.  (5c),  (9b),  and  (10)  into  Eq.  (1),  we  obtadn  for  the  linearized  current. 


j;(t)  =  ^  j;;  £  5^  exp[-.K(z?  +  -  aj))  \{C\  +  C75<)e-"‘ 

r  /  (*>  -  kv%  iKk(v°  1* 

-  /dwS*iV,At(w)exp[.'fc(a°  +  v^)  -  M(  -klfh^ 

-  ^  /  dw  u;A^(u;)  exp(ifc2y  —  iwt), 

Bo 


) 


(11) 


where  the  even  symmetry  of  S,  5*  =  5-fc,  has  been  used.  Here  constants  have  been 
absorbed  as  needed  into  the  constants  of  integration  C\  and  Cj-  Ignoring  aliases  again  and 
using  the  fact  that  the  initial  value  terms  (C}  -f- eventually  be  dominated  by 
terms  proportion2il  to  c"‘“*,  Im(a;)  >  0,  we  obtain  from  Ampere’s  Law 


j  du'  D+{k\^’]At{u}']e-'‘^'‘  = 

{§i)  dwB^(fc,w,A:\v°)A+(u;)e-‘<*'-‘)*^xp(-«(u;-(fc-fc>°^ 


Here, 


£?+(ifc',u;')  = 

k'W 


w'  -  fc'v? 


fc'v.ii  -  n  ’  (w'  -  fc'i;,||  -  n)2 


+ 


)• 


(13a) 
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is  the  linear  dispersion  funct' jn  that  would  result  if  each  of  the  particles  were  spread 
uniformly  in  physical  space,  and 

03. 

is  the  linear  coupling  coefficient  for  coupling  from  wavenumber  k  and  frequency  u;  to 
wavenumber  k'  due  to  scattering  off  particle  group  i.  Note  that  the  complex  exponen¬ 
tial  in  Eq.  (12)  forces  the  scattered  frequency  to  he  uj'  =  tj  —  (k  - 

The  coupling  process  may  be  diagnosed  by  formally  expanding  Eq.  (12)  in  orders  of 
(.V/.V,)"’.  (The  actual  expansion  parameter  is  not  in  general  (iV/JV,)“'  but  depends  on 
the  statistics  of  how  the  particles  were  loaded.  For  example,  if  the  particles  were  randomly 
loaded,  the  expansion  parameter  may  be  expectefl  to  be  of  the  form  a  constant  times 
(.V/A’j)“‘''^.  The  constant  will  usually  be  appreciably  greater  than  one,  since  typically 
the  resonant  particles  dominate  the  sum  over  i  in  Eq.  (12)  and  thus  effectively  prevent  the 
larger  number  of  non-resonant  particles  from  contributing  to  the  statistics.) 

The  effect  is  most  clearly  illustrated  by  considering  the  coupling  from  a  single  “pump" 
wave  with  designated  wavenumber  and  frequency  ko  and  uo-  The  only  equation  that 
survives  to  zero-order  in  (JV/JVs)~'  is  the  ko  component  of  Eq.  (12): 


D''’(ko,u;o)Ako(<*^o)  =  0 


(14a) 


which  just  indicates  the  pump  wave  must  satisfy  the  uniform,  linear  dispersion  relation  to 
zero-order.  The  k'  ^  k  equations  (coupling  equations)  first  appear  at  first  order: 


(Ao,  wo;  v°)>l+Ju-o)c-<'''-‘‘>)*-‘’ 


exp|-»(wo  -  [ko  -  fc')v'||)t).  (14b) 


There  is  thus  one  scattered  wave  for  each  group  of  pairticles  having  amplitude 


and  scattered  frequency 


w'*,  =  Wo  -  [ko  ~  k')v,ij. 


(15) 

(16) 


Equations  (14b)  and  (15)  and  definitional  Eqs.  (13a),  (13b),  and  (16)  serve  to  charac¬ 
terize  the  linear  mode  coupling  process  in  one-dimensional  AIC  simulations. 

Observe  first  the  complex  exponential  factor  involving  the  unperturbed  particle  group 
positions  z°  in  Eq.  (14b).  When  the  particle  group  positions  may  be  considered  to  obey 
random  statistics,  the  presence  of  this  factor  causes  the  right-hand  side  of  Eq.  (14b)  to 
scale  as  the  square  root  of  the  number  of  groups  ((iV/iV,)*^^)  as  suggested  above.  In 
contrast,  suppose  that  each  loaded  unperturbed  velocity  v°  is  represented  by  particle  groups 
i,, - ip, ... ,  ip,  i.e.,  v°  =  . . .  =  =  v°  which  are  uniformly  loaded  as  =  LzpjP  + 


Pump  waves  with 
their  radiation 


patterns 


Daughter  wave 


Loci  of 

I  0 

(j-kVi| 
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FIG.  1.  Linear  scattering  characteristics  of  AIC  simulations  in  {Re{uj),k) 
space.  Scattering  occurs  along  fan  lines  emanating  from  the  point  (0,0).  The 
slopes  of  the  fan  lines  correspond  to  the  parallel  velocities  present  in  the  ion 
distribution  function.  Radiation  patterns  are  shown  schematically  for  the  most 
effective  pump  waves,  i.e.,  those  with  the  largest  growth  rates  Im(u;).  An  exzunple 
of  linear  scattering  is  depicted  by  the  heavy  arrow.  Note  that  (1)  scattering  is 
from  a  mode  for  which  %  0  and  Im(u;)  is  large,  (2)  scattering  is  along  a 

fan  line,  and  (3)  (w',  k')  of  the  scattered  wave  lies  near  the  locus  of  k)  =  0. 


const,  in  the  periodic  simulation  system  (length  L,).  Then  for  any  pump  wave  (u'o,Ao) 
and  any  other  wavenumber  k',  the  sum  in  Eq.  (14b)  may  be  broken  up  into  sub-sums  each 
involving  a  single  value  of  v°  which  then  reduces  to  the  simple  sum  exp[— i(fc'  —  k)z°  ]. 

This  sum  always  vanishes  when  enough  particle  groups  are  loaded  per  loaded  velocity 
v°  so  that  |fc'  -  ko\Lx  <  2irP  for  the  maximum  wavenumber  difference  |fc'  -  fcol  allowed 
by  the  simulation.  When  the  number  of  groups  is  not  large  enough,  mode-coupling  has 
been  observed  for  wavenumbers  for  which  the  inequality  is  not  satisfied.  Details  of  this 
observation  will  appear  in  future  reports. 

We  next  turn  to  the  coupling  coefficients  D*(kQ,(j}o,k',y°)/D'^{k',ui')  as  they  appear 


in  Eq.  (15).  Clearly  the  coupling  will  be  most  important  between  those  modes  and  wave¬ 
lengths  for  which  \D'^(ko,ijJo,k' ,v^)\  is  large  while  |D‘''{fc',  w'^)!  is  small.  From  Eq.  (13b), 
we  observe  that  the  former  condition  is  most  easily  achieved  for  pump  waves  which  can  res¬ 
onate  with  some  portion  of  the  ion  distribution  function  (i.e.,  uq  —  —  fl  «  0),  while  the 

latter  condition  simply  requires  that  the  scattered  wave  be  of  wavenumber  and  frequency 
close  to  an  eigenmode  of  the  corresponding  uniform  system  described  by  D'*'(fc',w')  «  0. 
Figure  1  illustrates  the  consequences  of  these  two  conditions.  The  fan  of  lines  emanating 
from  {Q,k  =  0)  represent  the  locus  of  points  in  (Re(u;),  A)  space  which  are  in  resonance 
with  some  particle  group  in  the  ion  distribution  function.  There  is  one  line  for  each  value 
of  v°||  present  in  the  distribution.  The  fan  spreads  out  as  the  parallel  ion  beta  parameter 
is  increased,  with  the  line  representing  the  parallel  thermal  velocity  (for  a  bimaxwellian 
distribution)  intersecting  k  =  Cl/vy^  when  /Sqi  «  1.  For  /?q|  <  1,  the  uniform  system  disper¬ 
sion  curve  D'^{ko,ujo)  —  0  enters  the  thick  of  the  fan  of  lines  when  fco  >  ^/va-  Thus  the 
most  effective  pump  waves  for  linear  mode  coupling  may  be  expected  to  be  those  modes 
which  are  both  growing  with  significant  growth  rate  (since  the  coupling  also  depends  on 
the  amplitude  of  the  wave  through  the  factor  e*™l“)*)  and  have  ko  >  VI/va-  The  finite 
glm(u))t  gerves  to  smear  out  the  resonance  lines  in  the  fan.  The  dispersion  curve  may 
be  thought  to  lie  above  the  (Re(w),  k)  plane  for  Im(w)  >  0  or  below  it  for  Im(w)  <  0,  where 
the  third  dimension  corresponds  to  Im(cj).  The  coupled  modes  that  these  pump  waves  will 
drive  most  strongly  may  be  determined  from  the  wave-particle  scattering  relation  Eq.  (16). 
Since  the  line  connecting  (0,0)  and  pump  wave  (wq,  ko)  has  slope  which  also  appears 
in  Eq.  ( 16),  coupling  from  the  corresponding  pump  wave  would  expected  to  be  strongest  to 
modes  along  this  line.  More  precisely,  since  the  effective  pump  waves  typically  have  finite 
Im(u;o),  there  will  be  strong  coupling  to  any  available  modes  in  some  spread  about  the  line. 
Each  effective  pump  wave  may  be  thought  to  be  a  transmitter  in  (Re(u;),  k)  space  trans¬ 
mitting  with  frequency  Re(u;)  all  over  (w,  k)-land  with  some  radiation  pattern  determined 
by  the  details  of  Eq.  (13b)  and  the  magnitude  of  Im(wo).  The  radiation  patterns  of  the 
effective  pump  waves  are  shown  schmaticaliy  in  Fig.  1. 

The  best  “receivers”  (i.e.,  (w',  k')  of  the  scattered  wave)  clearly  lie  to  the  higb-k  side 
of  each  pump  wave,  as  illustrated  in  Fig.  1.  Here  the  receivers  are  more  in  line  with  the 
peak  of  the  pump  waves’  radiation  peaks.  Also  note  that  the  higher-k  pump  waves  have 
radiation  patterns  which  tend  to  line  up  more  closely  with  the  slope  of  the  dispersion 
relation,  implying  that  the  more  k  is  increased  the  tighter  the  coupling  becomes.  This 
effect  is  augmented  by  the  presence  of  the  factor  kk'  in  the  last  term  in  Eq.  (13b). 

We  conclude  then  that  linear  mode  coupling  due  to  the  discreteness  of  simulation 
particles  primarily  concerns  those  wavenumbers  for  which  roughly  k  >  O/uy^.  The  primary 
pump  waves  are  those  which  both  satisfy  this  condition  and  have  significant  growth  rates. 
The  most  strongly  scattered  waves  reside  to  the  higb-k  side  of  the  pump  waves  which  spawn 
them.  The  higher  k  modes  are  tightly  coupled  and  become  more  and  more  tightly  coupled 
as  k  is  increased.  In  stark  contrast,  the  lower  k  modes  (k  <  O/v^)  are  outside  both  the 
fan  of  Doppler-shifted  cyclotron  resonances  and  the  radiation  patterns  of  the  active  pump 
waves  and  so  do  not  participate  in  linear  mode  coupling.  All  of  these  conclusions  are 
consistent  with  simulation  results  obtained  so  far. 


E.  A  Model  of  the  Pla8ma>Sheath  Region  Including  Ion  Reflection 
Lou  Ann  Schwager  (Prof.  C.  K.  Birdsall) 

I.  Introduction 

Ion  reflection  at  the  collector  plate  affects  the  electrostatic  potential  and  transport  of 
the  incident  plasma.  Present  simulation  studies  indicate  that  increasing  ion  reflection  in¬ 
creases  the  electrostatic  potential  drop  across  the  sheath  region.  We  also  derive  theories  for 
the  dependence  of  potential  on  reflection  coefficient  R  for  cold  and  warm  ions.  Simulation 
results  compare  well  with  theory  for  small  R  at  the  ion  to  electron  mass  ratio  of  40.  For 
large  R  the  magnitude  of  potential  drop  from  the  computer  simulation  increases  beyond 
the  theoretical  value.  Also  the  oscillation  amplitude  widens.  As  we  observe  in  the  phase 
space  plots,  ion-ion  two-streaming  causes  an  instability  which  generates  large  amplitude 
oscillations  in  potential.  Hence  we  also  derive  a  dispersion  relation  for  this  two-streaming. 

II.  Simulation  Model 

The  computer  code  PDWl  modeb  this  interaction.  The  system  model  is  shown  in 
Fig.  1.  This  differs  from  the  model  in  the  previous  report  (QPR  I, II  1984).  Presently 
the  plasma  source,  at  the  left  boundary,  is  at  the  reference  potential.  The  plate  or  wall, 
at  the  right  boundary,  floats  to  the  wall  potential  The  plate  current  sums  to  zero. 
Previously  the  plate  was  at  the  reference  potential.  Also  presently  the  velocity  distribution 
of  the  reflected  ions  at  the  plate  matches  that  of  the  incident  ions  (specular  reflection). 
Earlier  the  user  specified  the  value  of  the  constant,  “reflected”  ion  flux.  The  present  model 
causes  all  reflected  ions  to  return  to  the  plasma  because  they  have  sufficient  energy  from 
the  incident,  accelerated  ions  to  overcome  the  sheath  potential.  Previously  with  a  half- 
Max”  ellian,  “re-emitted”  flux,  some  ions  could  not  overcme  the  barrier  and  so  returned 
to  the  wzdl.  Electrons  reaching  the  plate  are  absorbed.  The  phase  space  plots  in  both 


cases  clearly  indicate  this  effect.  As  before,  ions  and  electrons  returning  to  the  source  are 
thermalized  and  re-injected  (refluxed)  with  the  input  flux.  We  have  chosen  a  mass  ratio  of 
40  for  i?  =  0,  22,  44,  67,  and  89%.  * 


As  indicated  in  the  previous  report,  for  a  hydrogen  plasma  with  temperatures  of  30  - 
100  cV^  incident  on  molybdenum,  particle  reflection  coefiicients  vary  from  60%  to  80%. 
Reflected  ions  eiccount  for  only  10-15%  of  the  reflected  flux.  Neutrals  fill  in  the  remainder.* 
The  velocity  of  reflected  particles  depends  upon  the  relative  mass  ratio  of  the  incident  ions 
and  target  atoms.  The  ratio  of  reflected  to  incident  velocities  typically  is  about  90%.  We 
chose  to  simulate  the  wide  range  of  R  to  fully  understand  the  effect  of  ion  reflection  on 
potential. 

III.  Time-Independent  Analysis 
A.  Theory 


The  theory  for  the  dependence  of  potential  drop  on  ion  reflection  for  cold  ions  is  as 
follows.  Consider  an  infinite  wall  with  floating  potential  (an  open  external  circuit)  at  z=0 
in  contact  with  a  plasma  at  z>0.  Poisson’s  equation  describes  the  electrostatic  potential: 


with 


<P4> 

dx^ 


— (n*  -  n<  -  rir) 

«o 


(I) 


^(oo)  =  0 


where  rie  and  tii  are  the  densities  of  plasma  electrons  and  ions  and  nr  is  the  density  of 
reflected  ions.  Charge  neutrality  occurs  as  z  — » oo,  thus 

n,(oo)  =  n,(oo)  +  n,(oo)  =  n«,.  (2) 

*  Originally  R=20,  40,  60,  and  80%  was  input  but  a  counting  error  generated  the  actual 


numbers  given  above. 


Plasma  electroas  have  a  Maxwellian  velocity  distribution  with  temperature  T.  Their 
density  is 

ne(i)  =  f»oo  exp[c<^(i)/fcT).  (3) 

The  plasma  ions  are  cold  and  arrive  at  the  sheath  edge  with  an  energy  of  E  = 

Because  all  these  ions  entering  the  sheath  fall  freely  to  the  wall,  continuity  of  particles  and 
conservation  of  energy  determines  their  density; 

ni(a:)  =  -  ei^(z))~*/^,  (4) 

The  current  of  reflected  ions  depends  on  the  incident  current  at  the  wall  by  the  coef¬ 
ficient, 

R  =  —Jro/JiO'  (b) 

We  assume  these  reflected  ions  leave  with  the  incident  velocity  of 

VrO  =  ~ViO.  (6) 

No  net  current  flows  to  the  wall,  thus 

JtO  +  Jio  +  JrO  —  0-  (7) 

From  the  conservation  of  energy  for  reflected  ions  and  Eq.  (6),  we  have  that 

v’  =  2/M(E  -  (8) 

The  reflected  flux  is  constant  in  x.  With  this  and  Eqs.  (5)  and  (7),  the  reflected  ion  density 
becomes 
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We  then  substitute  Eq.  (8)  into  Eq.  (9)  to  determine  that 


Because  only  the  fastest  electrons  leaving  the  plasma  source  arrive  at  the  wall,  the  electron 


density  current  becomes 


JfQ  =  — erij 


OO 

i(m/2irkT)^^^  j  v  exp{—mv^ /2kT)dv 


where  u^mtn  =  —2e(l>o/m,  by  conservation  of  energy.  Integrating  we  find  that 

-Jto/e  =  fioo  exp(tlio)(kT /2irmy^^  (12) 

where  ^  =  ett>fkT.  We  substitute  this  into  Eq.  (10)  to  solve  for  reflected  ion  density, 

n,  =  n«exp(V>o)(^-^j  (£  -  e^)  ^  (13) 

We  find  the  incident  ion  density  in  combining  Eqs.  (2)  and  (4).  Next  we  substitute  the 
value  for  riroo  which  is  found  in  letting  <^—*0  in  Eq.  (13).  Thus 

The  three  densities,  in  Eqs.  (3),  (13),  and  (14),  are  then  placed  in  Poisson’s  equation 

(Eq-  (D), 

-^^=exp(tl-)-E‘/»(E-e.^)-‘/^  (15) 

erioo  dx^ 

At  the  plasma  source  where  x  -+  oo,  then  ^  » 0  and  —ddffdx  — ►  0  because  of  charge 
neutrality.  We  then  investigate  the  behavior  of  Eq.  (15)  as  z— *0.  Multiplying  by  dd>/dz 
and  integrating  from  oo  to  z  gives 


•  .  •  V  -  .  -  i.-*  .  ■  .  -  . 


fr: 


.'  /  /  / 


Expanding  e((>/kT  and  e^/ E  in  a  Taylor  series  and  retaining  the  6rst  three  terms  determines 
that 


kT  nao\dx) 


1  rl)e<l> 
4“F‘ 


(17) 


As  I— *00  then  the  electric  field  is  always  positive  smd  approaches  zero.  Hence  {d<(>/dx)^  >  0 
which  when  applied  to  Eq.  (17)  implies  that  as  i— »oo 


E  >  kT/2. 


(18) 


This  is  also  the  Bohm  condition  for  an  absorbing  wall. 

We  next  find  the  dependence  of  potential  on  R.  Re-arranging  the  zero  wall  current 
condition  and  using  Eq.  (12)  for  electron  current  density  at  the  wall,  gives 

/  jfcT  \ 

(1  -  R)Jio  =  efioo  e3cp(t&o)(  2^  j  •  (19) 


By  continuity,  ion  current  into  the  wall  equals  that  leaving  the  plasma  source.  We  find  this 
current  density  J^o  in  using  =  2E/M  and  Eq.  (14)  with  (f>—*0.  This  current  density  is 
substituted  into  Eq.  (19)  which  we  manipulate  to  obtain 


Using  E  >  kT /2  derives 

I 


(20) 


(21) 


The  efiPect  of  reflected  ions  on  the  heat  insulating  ability  of  the  sheath  is  evaluated  by 
calculating  the  net  energy  flux  Q  at  the  wall.  Each  incident  electron,  which  has  a  speed 
of  (r)  =  {2kT/irmy^^,  carries  an  energy  of  2kT  to  the  wall  by  thermal  efi'usion.  (See 


Appendix  for  derivation  of  thermal  effusion.]  The  incident  electron  density  is  one-half  of 


n,o-  An  incident  ion  brings  in  the  same  energy  as  a  reflected  one  removes.  This  energy  is 
(E  —  e<i)o}.  The  kinetic  energy  flux  to  the  wall  becomes 

Q  =  ^'^eo(v)(2kT)  +  riioVtoiE  -  e<i>o)  -  riroVroiE  -  e^io) 
or 

Q=-n^  exp(V-o)  —  (2kT  +  E-  e4>o).  (22) 

£  \  Km  j 

The  electron  free-flow  energy  flux  is 


Thus 

Q  =  Qsre.  *  G(R)  (24) 

where 

G(R)  =  ^  €xp(^o)  (2  +  -  V-o  j  • 

We  use  our  previous  results  from  Eqs.  (18)  (E  =  kT/2)  and  (21)  to  see  the  dependence  of 
G  on  R.  For  a  hydrogen  plasma  interacting  with  a  purely  absorbing  wall  then  G(0)  =  0.156. 
For  a  low  temperature  plasma  with  a  metal  collector  plate  Rmax^0.l2  (Reflected  ions  are 
15%  of  the  reflected  particles  which  are  80%  of  the  incident  ions.)  Then  G(0.12)  =  0.128. 
Hence  ion  reflection  improves  the  thermal  insulating  effect  of  the  sheath  at  most  by  18%. 
If  some  material  were  found  to  have  an  ion  reflection  coefficient  of  50%,  then  the  thermal 
insulating  effect  would  be  improved  by  60%. 

We  also  derive  the  dependence  of  potential  drop  on  reflection  coefficient  when  the  ions 
are  warm  and  the  ion  and  electron  temperatures  are  equal.  The  model  follows  the  previous 
derivation  for  cold  ions.  The  electrons,  primary  ions,  and  reflected  ions  have  a  Maxwellian 
velocity  distribution  with  temperature  T  at  the  plasma  source.  Here  charge  neutrality  also 


applies  (Eq.  (2)).  At  the  wall  the  reflected  ion  current  depends  on  R  and  the  incident 
velocity  (Eqs.  (5)  and  (6)).  The  wall  current  sums  to  zero  (Eq.  (7)).  The  electron  wall 
current  density  is  specified  with  Eq.  (12).  We  assume  that  all  ions  emitted  at  the  source 
reach  the  wall.  Similarly  all  ions  reflected  from  the  wall  arrive  at  the  source.  Hence 


V 

'\2irM  J 


*JrO  ”  CTlf 


/  kT 

'V25rM/ 


Substituting  Eq.  (25)  into  the  current  balance  of  Eq.  (19)  yields 

fioo  exp(V)o)  =  (1  -  •R)nioo- 


By  the  definition  of  R,  charge  neutrality,  and  Eqs.  (25)  and  (26)  then 


Woo  —  fliioo(l  "b  ^)- 


The  dependence  of  (t>o  on  R  is  found  when  we  substitute  Eq.  (28)  into  Eq.  (27): 


e<t>Q  _  1  1  +  i?  /  M  \  ' 

~  l-R\m) 


'  •  t.  "  u 

v.'.’.V.' 


B.  Results 


We  graph  both  results  (Eqs.  (21)  and  (29))  and  compare  these  with  simulation  in 
Fig.  2.  The  simulation  uses  a  plasma  source  of  Maxwellian  ions  and  electrons  at  equal 
temperatures.  The  source  sheath  accelerates  ions  above  the  sound  speed.  Thus  after  this 
acceleration  into  the  central  region,  where  charge  neutrality  results,  the  ion  temperature 
is  1.5  to  2  times  less  than  the  electron  temperature.  This  was  also  observed  in  the  phase 
space  plots  of  the  previous  QPR  for  a  purely,  absorbing  wall.  Thus  simulation  results  for 
4>w  fall  between  the  two  theories  but  much  closer  to  the  cold  ion  derivation. 
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Note  also  that  beyond  22%  the  amplitude  of  oscillation  increases  considerably.  The 
history  of  wall  potential  for  R  =  22,  44,  67,  and  89%  is  shown  in  Figs.  3  and  4.  (The 
previous  QPR  shows  potential  history  at  fi  ^  0%.)  As  R  increases,  the  onset  of  large 
amplitude  oscillations  occurs  earlier  in  time.  Also  the  amplitude  relative  to  time-averaged 
potential  increases.  The  reason  for  this  is  seen  in  the  phase  space  plots. 

For  R  =  22%  at  a  time  of  125,  phase  space  plots  for  ions  and  electrons  and  potential 
profile  are  in  Fig.  5.  The  electrons  have  a  quiet,  cutoff  Maxwellian  distribution.  The  ions 
stream  with  a  thermal  spread  and  follow  the  potential  profile.  Note  that  a  few  ions  have 
fallen  below  the  sharp  edge  velocity.  This  indicates  a  time-dependent  interaction  of  the  ion 
streams.  At  a  later  time  the  beam  velocity  of  each  stream  has  fallen  to  the  ion  sound  speed, 
defined  as  vs  =Ute(m,/me)i.  This  phase  space  configuration  shown  has  been  constant  for 
the  last  20%  of  the  simulation. 

For  R  =  89%,  the  phase  space  and  potential  plots  at  times  of  25,  37.5,  and  100  in 
Figs.  6,  7,  and  8  indicate  a  strong  interaction.  In  Fig.  6  the  ion  streams  begin  to  pinch 
together.,  presumably  in  some  form  of  two-stream  interaction.  Later  in  Fig.  7  vortices  have 
evolved.  Finally  in  Fig.  8  the  streams  have  collapsed  and  have  reached  a  new  equilibrium. 

In  each  case  the  electrons  appear  to  play  no  role  in  the  interaction.  Also  at  a  time  of 
100  the  number  of  particles  is  still  increasing. 

IV.  Time-Dependent  Analysis 
A.  Theory 

We  gain  a  better  understanding  of  the  ion-ion  two-streaming  in  solving  its  dispersion 
relation.  Assume  that  two  opposing  ion  streams  exist  far  from  a  wall.  The  primary  stream 
has  a  current  density  of  erxipVQ.  The  reflected  stream  has  a  current  density  of  -en^rVo.  In 
one  dimension  with  the  previous  deflnition  for  reflection  coefficient  R  and  continuity  then 

R  =  (30) 


-3 


By  charge  neutrality  the  total  ion  density  equates  with  the  electron  density  no  so  that 

nip  =  no/{l  +  R)  (31) 

"ir  =  no/[^/(l  +  fi)].  (32) 

Both  species  of  ions  are  cold.  Plasma  electrons  have  a  Maxwellian  velocity  distribution 
with  temperature  T. 

The  subsequent  derivation  follows  that  illustrated  in  Chen  for  the  electron  two-stream 
instability.^  We  look  for  electrostatic  waves  with  solution 


Ei= 


with  X  along  uq  and  k.  Next  linearize  the  equation  of  motion  for  the  ions  and  arrive  at  the 
perturbed  velocity  for  each  stream. 


-  ieE  (  1  \ 


a  =  p,r 


The  linearized  continuity  equation  for  the  ions  gives  the  perturbed  density 


riia  =  KnaVia/{ul  -  KVa) 


a  =  p,r 


The  Boltzmann  electrons  have  density 


where 


file  =  eno4>i/kT 


E  =  — IK^l. 


We  then  substitute  Eq.  (33)  into  (34)  and  Elqs.  (34),  (35),  and  (36)  into  the  linearized 
Poisson’s  equation, 

i/cE  = -^(n,, n,p  -  nie).  (37) 


The  resultant  dispersion  relation  is 


where 


and 


1  +  = 


2„2 


+ 


R 


1  4-  i?  [(w  —  KVq)^  (u)  +  KVoy 


Vg  =  kT /M 


Xp  =  (kT/q^no. 


(38) 


Note  that  we  recover  the  dispersion  relation  for  ion  acoustic  waves  when  tio  =  0  and  R  =  0. 

Consider  the  form  of  Eq.  (38)  as  iZ  =  0.  This  equation  becomes  quadratic  in  which 
gives  either  four  real  or  two  imaginary  and  two  real  solutions  for  each  real  k.  Thus  the 
value  of  K  for  marginal  stability  occurs  where  w  =  0.  Equation  (38)  then  yields,  for 

1  +  K^X])  « 

Thus  for  small  k,  such  that  kA^^Ii  then 


vs  >  Vo- 


B.  Results 

We  observe  the  change  In  instability  growth  as  vq  and  vs  change  their 
reiative  values  in  time.  Early  in  our  simulations,  vo  >  3v5.  Ion  velocity  then 
approaches  vs  within  the  time  required  for  a  thermal  ion  to  travel  to  and  from 
the  wail.  This  time  is  2£/vti  =  4/.158  =  25.  For  each  value  of  R,  the  ion  streams 
begin  to  pinch  when  vo^vs-  For  A  =  89%,  Figs.  9A  and  OB  show  the  histories 
of  electron  thermal  velocity  and  of  ion  beam  velocity.  Simulation  values  Just 
before  and  after  the  time  when  vo  ss  V5  are  substituted  into  Eq.  (38).  Figures 


lOA,  lOB,  and  IOC  illustrate  the  solution  with  these  values.  For  vq  >  vs  only  real 
roots  exist.  For  vq  <  vs  solutions  with  imaginary  u  emerge  indicating  growth  of 
the  instability.  Observe  the  pinching  of  the  streams  in  Fig.  8B  where  vq  <vs- 
Thus  Eq.  (36)  crudely  models  the  two-streaming  we  observe  by  predicting  the 
Initiation  of  the  instability.  However  the  growth  rate  predicted  {e.g.  at  89% 
the  maximum  value  of  Imaginary  u  is  about  1)  is  an  order  of  magnitude  greater 
than  that  observed.  This  is  because  the  dispersion  relation  assumes  cold  ions. 
Simulation  ions  have  a  thermal  spread  which  is  damping  the  instability. 

As  indicated  earlier,  small  values  of  R  provide  a  modest  two-stream  in¬ 
teraction.  From  the  uj-k  plots,  we  observe  that  as  i?  is  decreased,  the  value 
of  K  for  minimum  stability  with  vo  <  also  decreases.  In  our  simulations 
Hence  we  may  be  observing  a  modest  interaction  for  small 
R  because  our  system  length  is  too  short. 

Usually  the  initiation  of  the  instability  occurs  within  the  time  that  the 
slower  simulation  ions  have  travelled  across  2L.  Thus  instead  of  intializing  the 
system  empty,  we  may  try  to  begin  fully  loaded  with  Maxwellian  particles. 

V.  Suggestions 

Future  work  will  include  the  effects  of  system  length  and  simulation  ini¬ 
tiation  on  the  two-stream  interaction.  We  will  also  study  the  dependence  of 
potential  and  instability  growth  on  mass  ratio.  Diagnostics  for  temperature 
profile  and  particle  and  energy  flux  to  the  wall  will  also  be  implemented. 
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Appendix 


We  present  the  following  derivation  of  the  energy  transported  by  ther¬ 
mal  effusion  after  a  long  library  search  in  vain.^  Assume  that  particles  have  a 
Maxwellian  distribution  of  velocity  with  temperature  T  in  three  dimensions  z, 
y,  and  z.  We  are  interested  in  the  kinetic  energy  flux  along  z.^  This  kinetic 
energy  flux  is 


where 


and 


£.= 


Ill  ^rnv^Vifod^v 

—  oo-ooO 

oo  oo  oo 

/  /  J  fod^v 

—  oo-odO 


fo  —  A  exp  I 
rf  =  kT /m. 


'vl+vl  +  vl\ 
-2vf  )' 


=  n^27rv*j 


-3/2 


(Al) 


Thus 


OOOOOO  OOOOOO  QCOOOO 

^£x=  j  j  j ^mvlfod^v j  j  j ^mvlv^fod^v  +  j  j  j ^mvlvxfod^v{A2) 

— oo-ooO  — oo-ooO  -~oo-oo0 

“  ^xy  ^XM* 


Next  we  evaluate  each  term  of  Eq.  (A2), 


=  iitmv^AIn 


or 

=  {«)  kT  (A3) 


where 


=  2vmv^A/n 


or 

STy  =  \{v}kT.  (A4) 

Similarly 

£T.  =  \{v)kT.  (A5) 

Combining  Eqs.  (A3)- (AS),  we  find 

€^  =  {v)2kT.  (A6) 

Thus  the  average  particle  carries  an  energy  of  2kT  along  the  x  direction. 
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Figure  t.  Plasma  Model 
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Figure  5.  Profiles  in  (a)  electron  phase  space,  (b)  ion  phase  space,  and 
(c )  potential  for  22%  ion  reflection  at  equilibrium  ( t  *  125  ) 
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Figure  6.  Profiles  in  (o  )  electron  phase  space,  (b)  ion  phase  space,  and 
(c)  potential  for  89%  ion  reflection  at  t*  25 
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Figure  7.  Profiles  in  (a  )  electron  phase  space,  (b)  ion  phase  space,  ond 
(c)  potentiol  for  69%  ion  reflection  at  \-Z7.5 
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Figure  8  Profiles  in  (o )  electron  phose  space,  ( b)  ion  phase  spoce,  and 
(c)  potential  for  89%  ion  reflection  at  t«IOO 
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Figured.  Histories  of  (a) electron  thermal 
velocity  and  (b)  ion  beam  velocity  for  R*89% 
v»  *  40‘i  Vfh, .  Time  in  step  number  *  D  T. 
Velocities  are  averaged  over  the  central  half 
of  the  simulation  length. 
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F.  PLANAR  MAGNETRON  DISCHARGES 


Perry  Gray  (Prof.  C.K.  Birdsall) 
Amy  Wendt  (Prof  .  M.A.  Lieberman) 


This  is  a  new  project,  initiated  recently,  primarily 
by  Varian  Associates.  It  is  related  to  our  other  projects, 
through  the  study  of  the  general  problem  of  magnetized  sheaths. 


INTRODUCTION 

We  are  using  two  approaches  to  understand  planar  crossed  field  (mag- 
netron)  discharges:  experiments  and  computer  simulation.  We  are  proceeding 
with  the  design  and  construction  of  a  planar  discharge  device  wbdch  is  well 
instrumented  to  measure  the  space  and  time  dependence  of  potentials,  fields, 
densities,  and  velocities  as  a  function  of  a  variety  of  parameters.  We  are  also 
using  plasma  computer  simulation  for  a  similar  configuration  in  order  to  "meas¬ 
ure"  the  same  qualities.  This  interplay  between  experiment  and  simulation  is 
expected  to  be  as  productive  of  understanding  in  this  area  as  it  has  been  in 
other  areas  of  plasma  research  (e.g.,  fusion).  Part  A  of  this  report  describes  the 
experiment,  showing  the  design  and  giving  future  plans.  Part  B  describes  the 
computer  simulations  done  to  date  and  the  next  steps  in  modeling. 

The  past  year  (August  1983-August  1984)  is  our  initial  step  in  the  study  of 
crossed-field  discharges.  There  has  been  considerable  preparatory  effort,  exa¬ 
mining  the  physics  and  parameter  ranges,  designing  and  building  the  experi¬ 
ment,  and  using  our  existing  simulation  codes  on  rudimentary  models.  As  yet, 
there  are  no  experimental  results;  there  are  some  simulation  results  of  a  gen¬ 
eral  nature.  However,  results  directly  applicable  to  planar  magnetron 
discharges  are  yet  to  come. 

A  typical  magnetron  discharge  is  shown  in  fig.(l);  imagine 
this  model  rotated  about  a  vertical  axis  standing  on  its  left  side, 
so  as  to  obtain  the  so-called  racetrack  magnetron  discharge. 

A  EXPERIMENT 
1.  Summary  of  Progress 

Mechanical  design  of  the  experiment  has  been  completed,  and  machine 
shop  work  is  in  progress.  The  experimentad  chamber  and  end  plate  design  is 
shown  in  Figure!  •  *^0  8"  diameter  plates  are  mounted  on  the  axis  of  the  12" 
diameter,  30"  long,  aluminum  vacuum  vessel.  Each  plate  is  water  cooled,  insu¬ 
lated  from  the  vacuum  vessel,  and  czui  be  moved  along  the  axis  (by  means  of  a 
hand  crank  at  each  end)  while  the  system  is  under  vacuum  and  the  plasma  is 
present  Seven  2"  diameter  ports  in  the  plane  of  syrnmetry  provide  diagnostic 
access.  Axial  access  (through  each  8"  plate)  is  also  available.  Tlie  azimuthal  sym¬ 
metry  of  the  design  allows  comparison  of  two-dimensional  (r,z)  measurements 
with  two-dimensional  computer  simulations. 

Vacuum  design  of  the  experiment  has  been  completed,  and  all  new  equip¬ 
ment  necessary  has  been  purchased  and  received.  The  design  makes  use  of  an 
existing,  six-inch,  oil  diffusion  pump  system.  The  chamber  pressure  can  be 
varied  in  the  range  10“*  to  10  torr,  and  will  be  feedback  controlled  in  the  range 
10~*  to  1  torr  using  a  1-3/8”,  motor-driven,  butterfly  throttle  valve  and  a  capaci¬ 
tance  manometer  A  calibrated,  variable  leak  will  be  used  to  set  gas  flow  rates. 
The  first  experiments  will  be  done  in  argon  plasmas. 


The  cooling  and  electrical  system  designs  have  been  completed.  Most  com¬ 
ponents  (cooling  manifold,  etc.)  are  in  hand,  and  minimal  machine  shop  work 
will  be  necessary.  A  0-5  IcV,  0-3  A  dc  power  supply  for  the  discharge  is  in  hand  for 
the  initial  experiments.  A  higher  current  source  will  be  purchased  if  needed. 

Space  for  the  experiment  has  been  cleared  in  Room  111,  Cory  HalL 

Tlie  support 

structure  has  been  designed,  all  necessary  materials  have  been  ordered,  and 
most  have  been  received. 

The  electromagnet  design  has  been  completed,  and  is  shown  in  RgureJ.  All 
necessary  materials  have  been  ordered.  Machine  shop  work  will  commence  after 
machining  of  the  experiment  chamber.  The  coil  will  be  wound  with  104  turns  of 
0.125"  o.d.,  0.065”  i.d.  copper  tubing,  insulated  with  2  mil  mylar  tape.  The  coil 
resistance  is  estimated  to  be  0.057  0.  For  a  current  of  50  A  (5200  ampere-turns), 
the  ohmic  power  dissipation  will  be  140  W.  For  this  current,  the  magnetic  induc¬ 
tion  at  the  pole  faces  is  estimated  to  be  1.3  Tesla  The  cooling  water  flow  rate 
(two.  520  turn  sections  in  parallel)  has  been  measured  to  be  7.8  emVsec.  For  a 
magnet  current  of  50  A,  the  temperature  rise  will  be  small  4.4”  C.  The  magnet 
will  be  powered  using  an  existing  0-60  A.  0-12  V  power  supply.  At  50  A.  the  mag¬ 
net  voltage  will  be  2.9  V. 

Plasma  diagnostics  will  at  first  include  the  use  of  Langmuir  probes  and 
Faraday  cups.  Langmuir  probes  are  used  to  determine  ion  density,  electron 
temperature,  and  fioating  potential  within  the  plasma  Two  Langmuir  probes  and 
their  probe  driving  circuits  are  available  for  use  in  the  experiment.  Faraday 
cups  can  be  used  to  measure  ion  and  electron  fluxes  to  surfaces,  and  to  obtain 
information  on  the  energy  distributions  of  these  fluxes.  A  design  of  miniature 
Faraday  cups  for  the  cathode  and  anode  end  plates  has  been  completed. 

A  time-shared  link  to  the  existing  plasma  data  acquisition  system  is 
planned.  Nearly  all  the  equipment  required  is  presently  on  hand. 

2.  Futvn  f^cms 

Machine  shop  work  on  the  experimental  chamber  will  be  completed.  The 
support  structure  will  be  constructed,  and  the  chamber  and  vacuum  system  will 
be  assembled  and  tested.  A  dc  glow  discharge  will  be  established  and  used  for 
initial  checkout  of  machine  and  plasma  diagnostics. 

After  machining  of  the  chamber  is  completed,  the  electromagnet  will  be 
machined,  wound,  and  tested.  It  will  then  be  mounted  on  the  cathode  plate  and 
used  in  the  generation  of  planar  magnetron  discharges. 

Measurements  of  magnetic  field  B  will  be  made  using  a  Bell  615  Gaussme- 
ter,  which  is  available  on  short  term  loan.  Langmuir  probe  measiirements  of  ion 
density  n^,  floating  potential  if,  and  electron  temperature  T,  will  subsequently 
be  made.  A  Faraday  cup  will  be  mounted  on  the  cathode  plate  and  its  usefulness 
in  determining  the  local  ion  flux  and  energy  distribution  will  be  studied.  Small 
magnetic  field  (insulated  loops)  and  electric  field  (insulated  wires)  probes  will  be 
used  to  search  for  rf  fluctuation  activity  in  the  plasma.  Such  fluctuations  have 
important  implications  for  particle  and  energy  transport  in  the  discharge. 

The  time-shared  link  to  the  plasma  data  acquisition  system  will  be  installed 
and  tested.  'Hiis  link  will  make  it  possible  to  obtain  detailed  local  measurements 
of  B,  n^,  if,  T,,  and  rf  fluctuation  levels  (if  any)  as  functions  of  r  and  z  (the  sys¬ 
tem  is  azimuthally  symmetric).  If  Faraday  cup  measurements  are  successful, 
then  measurements  of  ion  and  electron  flux  as  a  function  of  r  can  also  be  made 
at  the  cathode  and  anode  plates. 
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This  second  phase  of  detailed  (r-z)  measurements  will  extend  into  the  third 
year  of  our  program  to  construct  a  well-instrumented  planar  magnetron 
discharge  and  to  study  its  properties.  The  objective  of  this  second  phase  is  to 
provide  the  data  on  particle  and  energy  densities  and  transport  coefficients 
required  to  compare  the  experimental  results  to  computer  simulation  results. 
Using  these  experimental  and  simulation  studies,  we  intend  to  construct  a 
model  of  the  planar  magnetron  discharge.  Initially,  a  point  model  (averaging 
over  r  and  z)  will  be  developed.  However,  the  goal  of  the  analytical  work  is  to 
develop  a  fully  two-dimensional  (r-z)  model  of  the  discharge,  as  a  function  of  the 
discharge  control  parameters  (pressure,  voltage,  power,  magnetic  field 
strength,  surface  materials,  and  geometry).  We  believe  such  a  model  will  be  very 
useful  in  understanding  how  the  discharge  scales  to  larger  sizes,  what  factors 
lead  to  non-uniformity  in  the  deposition  pattern,  and  which  pressures,  voltages, 
geometries,  etc.  maximize  the  deposition  rate. 

B  COHPUTEK  SmUlATlON 
1.  Summary  of  Progress 

We  have  begun  applying  computer  simulation  to  crossed-field  glow 
discharges  by  starting  with  a  simpler  model  which  is  better  understood.  How¬ 
ever,  new  physics  is  being  added  step  by  step.  The  current  model  is  planar,  one 
dimensional,  bounded,  and  magnetized.  Thousands  of  charged  particles  are  used 
to  model  the  electrons  and  ions,  all  followed  self-consistently  in  their  own  and 
the  applied  electric  and  magnetic  fields. 

The  state  of  the  art  in  many-particle  simulation  is  found  in  any  plasma  and 
computational  journals  and  in  our  new  book*.  The  simulations  are  initial-value, 
boundary-vedue  solutions,  including  sources  and  sinks,  complete  with  non- 
intrusive  diagnostics.  The  generally  presented  results  are  graphical,  both  during 
and  at  the  end  of  computer  runs  of  a  few  thousand  time  steps.  In  one  dimension, 
X.  Vg,  v^,  V,,  (which  is  favored,  as  being  most  economical  and  most  easily 
transferable  from  Berkeley  to  other  computer  systems)  moderate  runs  on  a 
CRAY-1  computer  take  5-10  minutes  (10,000  electrons  and  ions  for  10,000  time 
steps)  and  about  50  to  100  times  that  on  a  VAX  computer.  In  two  dimensions,  x, 
y.  Vf,  Uy,  1/,,  there  must  be  substantial  time  increases,  as  many  more  particles 
are  needed  to  fill  the  second  dimension.  There  is  no  lack  of  Id  and  2d  codes,  but 
all  must  be  modified  and  run  carefully  to  produce  believeable  physics. 

The  major  support  for  the  PTS  Group  is  for  fusion  plasma  research, 
from  DOE  and  ONR;  the  Varian  Associates  support  allows  us  to  complement  such, 
moving  into  plasma-wall  research  with  non-fusion  applications. 

The  most  thorough  study  of  a  bounded  plasma  done  by  us  this  past  year  was 
for  a  Q-machine  model;  one  end  plate  was  an  electron  and  ion  source  and  the 
other  was  an  absorbing  cold  plate,  which  could  be  biased  or  fioating.  Normally 
there  are  two  sheaths,  one  at  each  end,  with  a  zero  field  region  in  between.  One 
finding,  reversing  an  assumption  made  for  decades,  with  ion  rich  emission,  was 
that  there  was  no  electron  hole  (it  was  filled)  at  the  emitter  sheath;  this  was 
found  by  comparing  Vomm  and  from  simulation  with  predictions  from 

time-'mdependent  theory,  which  gave  different  answers  and  was  subsequently 


improved  to  have  the  "filled"  hole  (with  excellent  agreement).  See  Figure 
Another  finding,  quite  striking,  was  the  plasma  turbulence  seen  in  electron  and 
ion  phase  space  plots,  and  the  moving  double  layer  seen  in  potential  plots,  for 
moderate  positive  bias.  While  the  disruptive  currents  found  have  been  observed 
before,  experimentally  and  in  simulation,  it  is  believed  that  the  showing  of  all 
plots  simultaneously  (a  movie)  is  a  major  step  in  identifying  almost  all  of  the 
anomalous  behavior.  See  Figure  Simulation  results  provide  us  with  the  velo¬ 
city  distributions  and  transport  coeCTicients  in  the  plasma  and  sheath  regions  of 
the  discharge,  for  comparison  to  experimental  results. 
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A  step  in  the  direction  of  understanding  cros  sed-fleld  models  is  to  add  a 
static  magnetic  field  at  angle  y'  with  the  axis  (between  plates).  Our  initial  obser¬ 
vations  of  magnetized  plasmas  (and  also  by  Roland  Cbodura,  IPP,  Garching  bei 
Munchen,  FGR)  are  that  with  V*  <  90*.  the  usual  positive  plasma  is  observed 
(electrons  lost  to  the  walls);  however,  we  find  with  i/  =  90*  (jS  parallel  to  the 
plates),  the  plasma  becomes  negative  as  ion  gyroradii  exceed  those  of  electrons 
and  thus  ion  transport  exceeds  electron  transport  to  the  wall,  as  at  the  cathode 
surface  of  a  planar  magnetron  discharge.  For  ^  nearly  90*,  the  initially-negative 
plasma  soon  becomes  positive,  as  electrons  move  along  the  field  lines  to  the 
wall;  the  sensitivity  to  near  90®  needs  investigation.  Chodura's  theory  atnd 
simulations  indicate  a  magnetic  pre-sheath,  lengthening  as  V'"*90*;  he  has  no 
results  for  V'=90®.  We  plan  collaboration  with  Chodura  early  next  year,  partiadly 
to  resolve  the  "90®  sheath"  or  "side  wall  sheath,"  of  interest  to  both  fusion  and 
crossed-field  devices. 

Another  magnetized  model  is  shown  in  fig. 6.  A  ld3v  simulation 
has  been  started  in  the  xy  plane,  using  B  (x)  only,  and  we  are  looking 
at  a  slice  in  the  middle  of  the  linear  dipole  field. 


Futura  nms 

We  are  also  adding  binary  charged  particle  collisions,  scattering  with  neu¬ 
trals,  charge  exchange,  reflected  ions,  and  secondary  (and  other)  electrons  to 
the  code.  Each  of  these  steps  is  nontrivial.  We  expect  to  make  most  of  these 
additi'^ns  during  the  next  year. 

A  major  step  will  be  to  a  two-dimensional  model,  fully  crossed-field,  in  the 
cross-sectional  plane  of  the  typical  permanent  magnet  crossed-field  discharge 
used  in  sputtering.  This  model  will  be  developed  during  the  second  and  third 
years  of  our  studies. 
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Fig.  3.  Cross-section  of  cathode  configuration. 


Figure  5 •  Similar  to  Figure  4,  but  with  absorbing  plate 
biased  (by  a  battery)  with  a  moderate  voltage.  A  potential 
minimum  formed  at  earlier  time,  trapping  the  electron  vortex 
shown  In  electron  phase  space.  The  current  In  the  external 
circuit  shows  the  anomalous  effects,  repetitive,  with  both 
high  and  low  frequency  parts.  The  Ions  are  returned  to  the 
source  at  about  the  full  applied  potential  (many  times  v,.t,apnai 


L+Sl 


r. 


I 


G.  PARTICLE  SIMULATIONS  OF  THE  LOW-ALPHA  PIERCE  DIODE. 

Paper  presented  at  the  International  Conference  on  Plasma 
Physics,  June  27-July  3rd  1985,  Lausanne,  Switzerland. 


PAPTIULF  S  IML'LA  riO-.'S  rt<f;  LOW-A.’ •“'A  .'If'-.'L  U .  .-‘i 
T.L.  Ci'vstil  and  S. 

Plasn*  Theory  ar.d  S;muIatioh  Tr  •  it.  TKCS-EPL 
Vn  i  V.- r  a  i  tv  of  Ja  1 1  f.irn  la  ,  0«t<oU'V.  CA  '14'20 


Abat  fact .  Partule  s  i.t\u  lat  na  are  js'-d  stjo..-  me  evol-.it. or. 
of  small  initial  per tut bat  ion*  about  the  uniform  clasaioal 
P*erc€  diode  equilibrium  :n  the  paia.met«-r  ran9e  D  '  i ‘WpL.  v 
<  )•  .  Linear  yrorcrties  recovered  show  quan'itative  a  Teement 
With  rheoretital  predictions.  Lono-iime  nenavicr  is  (c  ;nd  to 
lepend  -n  a  and,  m  -.jne  case,  a*si  od  »he  ;r.:tiai  conditions. 

rudu-  t  I  -n .  Th*‘  ’rlassical'  Pierce  diode'  ^  is  charac- 
terired  by  e* ». e r -'a  1  ly  shotted  electrodes  lat  x«3.LJ,  a  -uni- 
ferm  bacx.jr  'ur.  1  jf  irrmoLil'.e  Ions  .  and  cold  olectfuns  injected 
constantly  jc  *•©  .lensJtv  and  vei,-ojty  v^ .  Out  of  the 

various  r.;u  i  1  ibr  1  jm  states  this  classical  p-.erce  diode  m.iy 
exhibit,^  the  jp.ifcrm  one  (3/hx«o,  3,'3t»0»  la  of  particular 
Importance  The  linear  e  luer  f  i  e<iuenc  le  6  asst.'Ciated  »»it.h  It. 
w  •  »  i*^  (with  u  ■  1.2,  ...I,  are  the  roots  of  tne 

characteristic  equation'  ' 

I  an  •  .  .  >  .  . 

e  Mn'»Usina  •  2incosQi  •  an  -  an'-2in  •  0, 

where  n«*ey«  a.nd  u  "w  L.'v  .  with  »»'n  e’  'i  ST.  The  potential 
P  po  pooe.^ 

pattern  of  eiqenmcde  „  .may  be  written  in  the  form'’ 

•  _^(x.tl  •  (xicoi  (u^^t»4l  •  ixl  sin  (U|^^t»4pexp(Yyt)  . 

All  linear  properties  are  conveniently  Classified  in  terns  of 
a.  The  a-values  considered  below  (O.b".  1.5*.  2.S^}  are  repre-* 
sentative  of  the  "archetypical*  a-intervals  (O.*),  {«,2al,  and 
i2*.3^>.  for  which  detailed  predictions  on  linear  behavior 
are  available. ^ ^  The  purpose  of  the  present  paper  is  to  check 
son*  of  these  predictions,  and  to  study,  in  the  linearly  -un¬ 
stable  cases,  the  nonlinear  evolution  as  well  as  the  "final" 
atates  of  the  classical  Piecca  diode. 


The  simulations  were  performed  with  the  new  bounded-plaecna 
code  PDWi^  tne  one'dlmenaional  electrostatic  mode.  Each  run 
was  initialiced  with  a  slightly  perturbed  electron  density 

profile  nis.t«0)  •  h^  ♦  ft  cce^**  {"poeltlve"  initiaiuetioni 
or  n{»,0»  •  n^  -  ft  cos-^  ("negetlve"  initiallxattonl .  where 
ft<<  n^.  In  all  cases,  the  number  of  simulation  particles  in¬ 
volved  wee  about  2000.  and  the  fields  were  celculated  on  a 
grid  of  12S  points  using  s  fmice-dlfference  Poisson  scheme. 

We  slso  choee  e/m^*i,  L-i.  and  as  useful  diagnos¬ 

tic  quantities  we  define  an  "internal  electrostatic  energy" 
and  Its  rate  of  change  by 

w*t)  -  ^j^dx^n^-n  ix,  t)  j  *  (X,  tl  and  >yltj  -  . 

respectively,  in  the  linear  regime,  we  only  consider  the  dom¬ 
inant  Iw-M  and  next-to-doninant  lw«2;  eigenmodes. 


.  Sw .  Theoretical  predictions:^  w^/w 
as  shown  In  fig.  'a 


-1.2271  jnon- 
/w  -5.412- 


oscillatory  decay) 

2.5461  (oscillatory  decay).  Results  of  e  sinulstion  run  with 
positive  iniciailzation  ere  shown  in  Figs.  ib-d.  From  Flge. 
lb  and  ic.  the  dominart-eigennode  stage  (in  which,  by  defini¬ 
tion.  all 
other  ei- 
genmodet 
are  negli¬ 
gible  rela¬ 
tive  to  the 
dominant 
one)  Is 
seer,  to  set 
4  in  a 


l«i 


^ract ically 
zero,  thus 
ronf  1  rising 


1  *  T  .••■-■ireticai  {  icdictluns;  '  *  0.l525i 

"ibcillatory  •jrow^r.,,  (not  shown  hrr-i;  •  2 .  3 1  iO  - 

O.o6n9i.  Thfse  .meur  ,vr&per  c  i  were  .-rrjf-ci  q  js:  1 1  tat  i 
as  described  tot  uRO.b".  T.he  global  evoluti-’h  of  tr.e  '  unstable ! 
system  f*)!  :-i-a<f<ve  ini  t  la  1 1  rat  ; is  i-'.-cu.-'-' r-,c«d  ty  Fig.  2>. 


according  to  which  the  diode  reaches,  at  t«)0.  a  new,  non-uni¬ 
form,  stable  d.c.  equilibciom  with  the  paiabclic  potehiial 
profile  sr.i.wr  in  rig.  2b.  Tr.is  is  in  sharp  contratt  to  th« 
case  of  initialisation,  in  wn.ch  the  final  state  Is 

characterized  by  large-amplitude  relaxatlo.’i  oecillatlone  (Fig. 
3a).  The  relatec  potential  profile  basically  aaintains  the 
shape  shown  in  Fig.  lb.  but  oecillatlone  of  the  virtue,  ceth- 
ode  near  the  injection  plane  periodically  interrupt  the  flow 
of  electrons  to  the  right,  thus  giving  rise  to  large-aapl itude 
OBclIlations  ih  the  external  current. 


4.  o  -  2. St.  Theoretical  predictions;^  w^^/wp  ■  0.1319  • 
0.07501  (oscillatory  growth),  ♦^^(x)  and  ♦^J(x)  (not  shown 
here);  •  1.9203  -  0.374di.  The  predicted  properties  of 

the  dominant  eigenmode  were  again  verified  as  before,  but  the 
next-to-dominant  elgenf requency  extrected  was  too  inaccurate 
to  be  conclusive.  Figure  4a  shows  the  global  evolution  for 
positive  initialisation,  which  again  Leads  to  a  state  of  sig¬ 
nificant  noise.  What  Is  new  is  that  there  are  now  tmo  oscil¬ 
lating  virtual  cathodes  (Fig.  4b),  which  alternatingly  reflect 
electrons  and  thus  give  rise  to  s  feirly  complex  distribution 
in  phsse  space  (Fig.  4c).  The  same  type  of  final  state  was  ob¬ 
served  for  negative  Initlailaation. 


*)Work  aupported  by  DOE  Crant  OE-AT03-76ET53604,  ONR  Grant 
NOOOI4-77-C-OS76,  and  Auatrian  Research  Funds  Grant  E-18/03. 
All  computations  were  performed  at  NKFECC,  Livermore. 

* ) Prrmannnf  .address;  Inetitute  for  Theoretical  Physics,  uni¬ 
versity  of  Innsbruck,  A-6020  Innsbruck.  Austria. 

'j.R.  Pierce,  J.  Appl.  Phys.  721  »1944). 

^O.B.  Godfrey,  Reports  AHPC-R-282  (1  981),  MRC-N-212  (1982). 
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Paper  presented  at  the  Second  Symposium  on  Plasma  Double  Layers  and 
Related  Topics,  July  5/6,  1984,  Innsbruck,  Austria. 


H  .  Theory  and  Simulation  of  Ion  Acoustic  Double  Layers 

If.  y.  f&m  and  T.  L.  Crystal 
E.R.L.,  University  ol  California,  Berkeley.CA.  94720 


1.  Simulations  of  Ion  Acoustic  Double  Layer 

All  simulations  are  done  with  a  Id  axially-bounded  electrostatic  PIC  code. 
In  all  of  our  simulations,  we  have  used  the  same  mass  ratio  (m^/m,  =  100);  the 
time  step  is  u,  it  =  0.2.  Initially,  the  simulation  plasma  density  is  uniform  in 
space.  The  ion  and  electron  distribution  functions  are  both  Maxwellian;  the 
ions  are  cold,  7)  «  7^,  and  the  electrons  are  drifting  relative  to  them  with  drift 
velocity  v^.  This  relative  drift  between  the  electrons  and  ions  constitutes  a 
current  and  can  result  in  instability  depending  on  v^. 

(a)  Simulations  of  current  driven  systems 

In  our  simulations  of  constant  current  driven  system,  we  have  found  that, 
weak  ion  acoustic  double  layers  can  be  formed  even  in  a  very  short  system 
(BOX,),  and  even  when  the  electron  drift  velocity  is  small  compared  to  previous 
simulations  (v^  =  O.ASvu, ’^4.5  6',);  the  double  layer  formation  mechanism  is  still 
based  essentially  on  ainpliflcalion  of  a  small  negative  potential  dip,  due  to 
reflection  of  electrons. 

Using  a  temperature  ratio  r  =  T,/  =  20  and  plasma  parameter 

nX,  s  100,  the  system  plasma  is  loaded  uniformly  in  space;  there  is  thus  no 
electric  fleld  initially.  From  Figure  1  '  note  that  by  time  u,  t  =  4B0  a 

small  negative  potential  dip  has  developed  that  is  associated  with  an  ion  phase 
space  distortion  as  well  as  with  an  ion  density  dip.  From  subsequent 
"snapshots"  of  this  potential  dip,  it  is  seen  to  be  moving  with  nearly  ion  acoustic 
velocity  (C,).  Therefore  it  can  start  to  trap  those  ions  which  are  resonant  with 
the  structure  (i.s. ,  ions  in  the  positive  velocity  tail  of  the  distribution),  and  an 
ion  hole  starts  to  form  there.  At  the  left  side  of  the  growing  potential  dip,  elec¬ 
trons  with  velocity  slightly  greater  than  the  ion  acoustic  velocity  contribute  to 
the  structure's  growth;  their  velocity  distribution  (in  the  moving  frame  of 
potential  dip)  has  positive  slope,  so  that  they  give  up  their  energy  to  this  poten¬ 
tial. 

As  the  negative  potential  structure  grows  and  decelerates  due  to  this  elec¬ 
tron  reflection,  the  dip  is  able  both  to  trap  more  ions  and,  at  the  same  time,  to 
reflect  more  electrons.  "Positive”  momentum  transfer  due  to  electron 
reflection  leads  to  "deceleration"  of  the  ion  hole,  and  the  potential  thus 
appears  to  have  a  negative  effective  mass  as  well  as  negative  effective  charge. 
This  deceleration  and  growth  of  the  potential  dip  can  lead  to  increased  ion 
trapping,  because  the  structure  sees  more  densely  populated  ion  distribution 
as  It  decelerates.  This  electron  reflection  causes  the  asymmetry  of  potential 
due  to  more  electron  density  buildup  at  the  left  hand  side  of  potential  dip  than 
that  of  right  band  side.  At  time  u,  (  =  BBO.  an  ion  acoustic  double  layer  is  well 
developed,  with  a  negative  net  potential  jump  |df)|  =  0.3  F,  over  a  distance  of 
about  lOX,. 


One  reason  for  initiation  with  the  low  drift  velocity  of  our  simulation  is  the 
condition  of  constant  current  injection  as  opposed  to  (tie  previous  simulations 
using  decaying  current  injection.  Of  cource,  there  is  a  similarity  between  our 
system  (with  constant  current  condition)  and  previous  long  periodic  simula¬ 
tions  with  decaying  current  conditions;  because  the  system  length  is  long  in 
the  periodic  simulation  and  because  of  periodicity  (a  particle  leaving  at  one 
boundary  is  replaced  at  the  opposite  boundary  by  an  incoming  particle  with  the 
same  initial  velocity),  early  in  time  there  is  noarly  constant  current  coming  in 
and  going  out.  This  nearly  con.stant  current  in  a  long  periodic  simulation  acts 
as  a  source  of  energy  which  leads  to  the  formation  of  weak  ion  acoustic  double 
layers  by  the  reflection  of  electron  current.  However,  our  bounded  simulation 
system,  with  constant  current  injection  has  more  energy  available  than  that  of 
a  long  periodic  simulation  system,  resulting  in  a  lower  threshold  drift  velocity 
than  that  of  a  periodic  system. 

In  all  of  our  simulations,  the  negative  potential  dip  always  moves  in  the 
same  direction  (t.e. ,  that  of  the  electron  drift  velocity).  This  is  because  right 
going  ion  acoustic  waves  are  preferentially  amplified  by  the  inverse  electron 
Landau  damping  but  left  going  ion  acoustic  waves  are  severely  attenuated  by 
the  Landau  damping.  The  early  growth  and  deceleration  of  the  potential  dip 
can  be  understood  qualitatively  as  a  simple  momentum  exchange  between  the 
structure  and  the  particles.  The  rate  of  change  of  ion  acoustic  double  layer 
momentum  should  be  approximately  equal  to  the  negative  of  the  rate  of  change 
of  particle  momentum;  this  can  be  expressed  by  considering  the  reflection  of 
electrons  and  ions  by  the  ion  acoustic  double  layer  as  follows(see  Fig.  2.(b)): 

jp  t'.i-J’.l*  t^.i)* 

- /  du  i/*/,(A/  +  u)  +  Zm,  f  dvv*f,(U-v) 

**  •4  % 

-  2ms  f  otw  u*  fi{M  +  u)  .  (1) 
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where  p,  =  *  f-Pn  -  •  Pi  =  — -yi  and  JU  represents  the  velocity  of  ion 
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acoustic  double  layer.  It  is  straightforward  to  write  down  the  small  amplitude 
limit('yi  ,  '^  «  1)  of  the  above  equation. 

The  first  term  in  Gq.(l)  represents  the  positive  momentum  lransfer(i.e. 
electron  momentum  loss)  due  to  electron  reflection  from  the  left  hand  side  of 
ion  acoustic  double  layer.  The  second  term  in  Eq.(l)  represents  the  negative 
momentum  transfer  due  to  electron  reflection  from  the  right  hand  side  of  ion 
acoustic  double  layer.  The  third  term  in  Eq.(l)  represents  the  positive  momen¬ 
tum  transfer  due  to  ioer  due  to  ion  reflection  from  the  left  hand  side  of  ion 
acoustic  double  layer. 

Initially  there  is  no  appreciable  asymmetry  in  the  potential  structure  (i.c. , 
ywO)  and  dPi^/  dt  >  0  since  electron  drift  is  in  the  right  direction  and  thus 

^  paoUivm  momentum  transfer(i.e.  electron  momentum 

loss)  by  electron  reflection,  the  ion  acoustic  double  layer  having  negative 
effective  mass  decelerates;  the  negative  effective  mass  results  from  the  ion  den¬ 
sity  dip  associated  with  the  ion  acoustic  double  layer.  As  potential  asymmetry 
develops  by  electron  reflection(i.e.  -y  >  0),  ion  reflection  starts  come  into  play 
and  the  second  term  starts  to  compete  with  the  flrst  term.  For  the  cold  ion 
case,  the  third  term  can  be  considered  negligible,  and  ion  acoustic  double  layer 
can  receive  net  positive  or  negative  momentum  transfer  and  thus  decelerate  or 
accelerate,  depending  on  the  relative  magnitude  of  fif  and  y,  the  velocity  of  the 


ion  acoustic  double  layer  and  the  electron  distribution  function  at  both  sides  of 
the  ion  acoustic  double  layer. 

(b)  Simulation  of  voltage  driven  systems 

For  the  constant  voltage  driven  system  with  current  injection  (ip  =  0.5), 
we  have  found  results  that  are  similar  to  the  above  current  driven  simulation. 
An  important  new  feature  however,  is  that  the  ion  acoustic  double  layer  now 
alxuaya  appears  Arst  near  the  left  side  of  wall  and  develops  more  quickly  than  in 
the  earlier  current  driven  simulation  (Fig.  ^  and  3).  In  addition,  our  constant 
voltage  driven  simulation  requires  a  lower  electron  drift  velocity  (0.2uot)  for  the 
formation  of  ion  acoustic  double  layers  than  was  necessary  in  the  constant 
current  driven  simulations. 

In  our  constant  voltage  driven  system  with  current  injection  at  the  boun¬ 
dary,  an  applied  potential  across  the  system  acts  to  increase  the  ef  fectvue 
drift  velocity  of  electrons.  Since  ions  injected  at  left  boundary  see  a  potential 
barrier  due  to  the  applied  potential,  they  will  be  reAected.  thereby  forming  an 
ion  phase  space  density  dip,  which  contributes  to  the  formation  of  a  negative 
potential  dip  near  the  left  side  of  the  wall(Fig.  3(a)}.  Therefore,  wc  expect  to 
And  formation  of  ion  acoustic  double  layers  at  a  lower  electron  drift  velocity 
and  near  the  left  side  of  the  wall. 


2.  Theory  of  loa  Acoustic  Double  Layers 

Using  our  graphic  method,  we  shall  show  that  there  is  a  'criticaT'  velocity 
for  the  existence  of  ion  acoustic  double  layers,  which  is  smaller  than  the  value 
reported  in  previous  papers,  and  that  there  arc  maximum  amplitudes  of  the 
potential  dip  corresponding  to  the  drift  velocities  exceeding  critical  velocity, 
ife  also  And  that  the  net  potential  jump  across  the  ion  acoustic  double  layers  is 
determined  by  the  temperature  difference  between  two  plasma  regions. 

From  the  simulations  we  see  that  the  electron  kinetics  are  important  and 
our  theory  accordingly  must  start  with  a  reasonable  yet  tractable  kinetic  elec¬ 
tron  modtl  which  uses  all  three  constants  of  motion.  To  this  end,  we  introduce 
the  following  "modiAed  Schamel"  type  of  electron  distribution  function,  having 
three  electron  components: 
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where  c  =  v*  -  2^.  The  Arst  component  is  the  "free”  (or  untrapped)  group  ol 
electrons.  They  make  up  the  bulk  of  the  electrons,  and  are  modeled  here  as  a 
drifting  Maxwellian  (unction  at  p  =  0;  as  given,  their  temperature  has  been  nor¬ 
malized  (t.e. .  to  T,/  •)■  ReAected  (or  "trapped**)  particles  populate  the  regions 
either  z  >  z,,  or  x  <  on  each  side  of  the  double  layer,  and  cannot  communi¬ 
cate  with  each  other;  therefore,  we  have  introduced  two  separate  temperatures 
(ft  and  6)  and  two  normalization  constants  (/g  and  / 4)  for  them.  Here,  electron 
velocity  (v),  the  potential  (p).  and  the  distance  (z)  are  normalized  to  the  elec¬ 
tron  thermal  velocity  (T./m,)^  ,  the  ''free**  electron  temperature  T,/t,  and 
the  electron  Debye  length  X,  =  (T./AimgC*)**,  respectively.  6  represents  simply 
the  Heaviside  step  function. 
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Since  we  are  interested  in  describing  ion  acoustic  double  layers,  we  shall 
use  a  fluid  formalism  to  describe  the  essentially  cold  ions: 

(n<)t  +  =  0  .  u,  u  u,  +  p,  =  0  (3) 

Considering  a  time  stationary  situation,  the  above  fluid  ion  model  yields 
TVi  u  =  no  Uo  and  uV  2  +  ^  =  £"0  where  no.  Uo  and  £"0  are  constants. 

From  these  two  species  models,  the  corresponding  densities  for  electrons 
and  ions  in  the  ion  acoustic  double  layer  frame  can  be  found  analytically  by 
velocity-integration  and  Poisson's  equation  may  be  written  by  introducing  a 
Sagdeev  potential  V(if)  as  follows: 

VxM  = 

The  Sagdeev  potential  is  found  from  the  obvious  first  integral  of  Poisson’s  equa¬ 
tion  with  appropriately  chosen  boundary  conditions. 

Numerically  solving  our  nonlinear  eigenvalue  equation,  we  have  been  able 
to  find  that,  in  fact,  there  is  a  “critical"  velocity  for  ion  acoustic  double  layer 
solutions  to  exist.  Arguing  that  the  smallest  visibl*  potential  dip  in  a  simulation 
would  be  of  order  the  corresponding  “critical"  velocity  for  this  ampli¬ 
tude  solution  to  exist  is  found  to  be  U.4'U|)i.  This  value  is  considerablly 

smaller  than  has  been  reported  in  other  papers,  and  much  closer  to  our  simula¬ 
tion  result  (i.e. .  that  i  was  necessary  before  the  current  driven  dou¬ 

ble  layer  would  form).  See  Fig.  4. 

Examining  solutions  to  our  nonlinear  eigenvalue  system,  we  have  found 
that  there  are  maximum  amplitude  limits  for  the  negative  potential  dipj-^i): 
these  depend  on  the  electron  drift  velocity  exceeding  the  critical  value  (Fig.  4.). 

We  have  also  calculated  ion  drift  velocities  in  the  frame  of  ion  acoustic 
double  layers  and  found  that  the  usual  Bohm  condition  is  not  met  in  the  case  of 
ion  acoustic  double  layers.  In  fact,  the  ion  drift  velocity  decreases  (below  C,) 
as  both  amplitudes  of  negative  dip  (-V'i)and  net  potential  drop  (+■#)  increase 
(see  Fig.  5  ).  in  confirmation,  the  simulation  DL  potential  minima 
are  always  observed  to  move  subsonically (Fig.  1  and  3). 

With  regard  to  the  common  identification  made  between  Ibe  ion 
acoustic  soliton  and  our  double  layer  solutions,  we  point  out  that  the  velocity  of 
the  usual  ("srefactive,  having  negative  polarity)  ion  acoustic  soliton  increases 
with  increasing  amplitude;  this  character  is  in  direct  conflict  with  our  earlier 
observation  that  the  ion  acoustic  double  layer  slows  down  as  it  grows.  It  is 
important  to  note  that  net  amplitude  of  an  ion  acoustic  double  layer  correlates 
directly  with  the  temperature  difference  between  the  two  plasmas:  the  net 
potential  drop  (V')increases  with  increased  temperature  of  reflected  electrons 
on  the  high  potential  side(see  Fig.2  (b)). 


We  would  like  to  thank  to  Prof.  C.  K.  Birdsall,  Dr.  M.  Hudson.  Dr.  S.  Kubn.  Dr. 
W.  Lutko,  Dr.  J-P.  Lynov,  Dr.  V.  Thomas.  We  are  grateful  to  Mr.  N.  Otani  and  Mr.  W. 
Lawson  for  their  computer  simulation  code.  This  work  was  supported  by  DOE 
Contract  DE-AT03-7aETS3064  at  ERL.  University  of  California.  Berkeley. 


SECTION  II:  CODE  DEVELOPMENT 


A.  Uniform  Number  Generation  (for  Quiet  Starts) 

N.  Otani 

We  describe  a  simple  pseudo-random  number  generator  based  on  n-  roman, 
reversing  .  To  each  integer  index  i,  a  real  funber  /  is  assigned  in  the  interval 
(0,1)  according  to  the  following  prescription: 


1)  Write  i  in  base-7i: 


i  =  2 

psQ 

Simply,  to  the  integer  i  expressed  in  base-n 

i  =  mp  ■  •  ■  rngmi-O  (basen) 


the  fraction 


/  =  O.mjTng  •  •  •  Tnp  (basen) 


is  assigned. 

For  the  purpose  of  loading  particles,  the  algorithm  implemented  need  not 
be  fast;  the  one  shown  here  simply  uses  the  above  definition: 


subroutine  nrev  (i,  fract,  n) 
c 

c  Returns  a  fraction  between  0  and  1  representing  the 
c  n -reversed  number  corresponding  to  i. 


*See  C.  K.  Birdsall  and  A.  B.  Langdon,  naama  nysies  via  Cbmputar  Sirr^lation,  McCraw  Hill, 
19B5,  Section  lS-5,  p.  393. 


!  i  index  to  be  n -reversed 

;  fract:  number  between  0  and  1  representing  the 
integer  in -reversed. 

;  n  :  the  integer  i  will  be  reversed  in  base  n. 
integer  i,  n 
real  fract 


r .  • 

[:■ 


J  =  1 

fract  =  0 
power  n  =  1. 

10  jnext  =  j/n 

powern  =  powern/n 
fract  =  fract  +  (j-inext*n)*powern 
if  (jnext  ,eq.  0)  return 
j  =  jnext 
go  to  10 
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(3)  N.  F.  Otani,  "Hybrid  simulations  of  the  gravitational  instability  in  the  pres¬ 
ence  of  a  large-amplitude  RF  wave." 
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SIMULATION  OF  OBLIQUELY  PROPAGATING 
BERNSTEIN  WAVES*. 

Wm.  S.  Lawson,  U.C.  Berkeley,  Electronics  Research  Laboratory. 

A  modification  of  the  code  PDWl  is  being  used  to  simulate  the  boun¬ 
dary  value  problem  for  Bernstein  waves  propagating  slightly  off- 
perpendicular  to  the  magnetic  field.  The  numerical  experiment  is 
modeled  after  the  Q-machine  experiment  of  Armstrong  et  al.ll].  In 
the  simulation,  Bernstein  waves  are  excited  at  one  end  of  the  simula¬ 
tion  (which  is  assumed  to  be  symmetric  about  the  exciter).  The  exciter 
is  a  set  of  three  grids,  the  outer  two  being  tied  together.  Imposing  a 
voltage  between  the  inner  and  outer  grids  at  a  particular  frequency  gen¬ 
erates  Bernstein  waves  which  are  seen  to  propagate  from  the  exciter 
across  the  system.  Wave  numbers  and  spatial  damping  rates  (due  to 
Landau  damping)  can  be  accurately  measured  and  compared  with  linear 
theory.  Effects  observed  include  the  onset  of  non-linearity  at  electric 
potential  amplitudes  of  only  a  few  percent  of  thermal  energies. 

ll]  Armstrong,  Rasmussen,  Stenzel,  and  Trulsen,  Physics  Letters, 
8SA(1981)281-284. 

•  Work  performed  under  Dept,  of  Energy  contract  DE-AT03- 
76ET53064  using  the  NMFECC  computers  at  Livermore. 


PLASMA  DIODE  SIMULATION:  THE  PDWl  CODE 
C.K.  BIRDSALL,  Wm.S.  LAWSON,  T.L.  CRYSTAL,  S.  KUHN.  N.F. 
OTANI,  A.B.  LANGDON.  E.R.L.,  Univ.  ofCaltf.,  Berkeley,  CA  94720.* 
PDWl  solves  for  electron  and  ion  motions  in  the  self-consistent  elec¬ 
trostatic  field  of  a  ld(l,2,3v)  plasma  device  having  external  series  RLC 
circuit  as  shown  below.  The  code  and  complete  documentation  are  now 
available.  Among  many  options,  one  may  model  open  circuit,  short  cir¬ 
cuit,  fixed  potential  drive,  constant  current  drive,  or  a  full  RLC  circuit. 
Standard  code  options  include  various  particle  injection  schemes  (e.g. 
drifting  cutoff  Maxwellians)  from  either  or  both  electrodes,  and  a  static 
magnetic  field  Bq  at  an  angle  to  the  system  axis.  Recent  applications, 
include  the  Pierce  diode,  weak  and  strong  double  layers  both  with  and 
without  Bq,  Bernstein  waves,  and  magnetic  sheaths  and  first  walls. 
PDWl  was  produced  for  and  during  a  Plasma  Device  Workshop. 


•  Work  was  supported  by  DOE  contract  DE-AT03-76ET53064  and 
ONR  contract  N00014-77-C-0578.  Computations  were  done  at 
the  NMFECC,  Livermore. 


Hybrid  Simulations  of  the  Gravitational  Instability  in  the 
Presence  of  a  Large* Amplitude  RF  Wave.  NIELS  F.  OTANI,  Elec¬ 
tronics  Research  Lab,  U.  of  Calif,  Berkeley,  CA  ^4720.’  A  2  l/2-(l 
particle-ion,  fluid-electron  nonlinear  Darwin  simulation  code  has  been 
developed  for  the  study  of  the  behavior  of  the  gravitational  instability 
in  an  equilibrium  modified  by  ponderomotive  and  other  effects  due  to 
the  presence  of  a  large-amplitude  RF  wave  in  the  ion-cyclotron  fre¬ 
quency  range.  No  change  is  observed  in  the  linear  growth  rate  of  the 
instability  when  the  RF  wave  is  loaded  as  a  long-wavelength  (ki—LH’'), 
compressional  Alfven  standing  wave  in  a  cold  plasma  with  initial  ampli¬ 
tude  e^Vl£lV4m^<u^g*-0(l);  however,  this  type  of  wave  lacks  the 
usual  ion-cyclotron  resonance.  A  theoretical  prediction  of  the  second- 
order  RF-wave-induced  d.c.  current  in  the  ki.x  Bo-direction  in  a  uniform 
plasma  agrees  with  results  obtained  from  simulation,  with  assessment 
of  the  current's  stabilization  properties  still  pending.  Further  simula¬ 
tions  employing  waves  exhibiting  the  ion-cyclotron  resonance  (i.e., 
finite- A(ii  Alfv'n  waves  or  ion-Bernstein  waves)  as  well  as  antenna- 
launched  waves  are  anticipated. 

‘Work  supported  by  DOE  Contract.  No.  DE-AT03-76ET53064;  compu¬ 
tation  facilities  provided  by  NMFECC. 


SIMULATION  OF  CHARGED  PARTICLE  REFLECTION  AND 
EMISSION  EFFECTS  IN  THE  PLASMA-SHEATH  REGION,* 

L.A.  SCHWAGER,  U.C.  BERKELEY,  E.R.L.,  CA  94720 

The  plasma  sheath  region  is  modelled  by  a  one- 
dlmenslonal,  time-dependent  system  using  PDWl,  which 
simulates  the  Injection  of  particle  Ions  and  elec¬ 
trons  from  a  plasma  region  Into  a  grounded  plate. 

At  this  plate  electrons  and  Ions  are  reflected, 
emitted,  and  absorbed.  Current  work  Indicates  that 
a  reflected  Ion  current  Increases  the  plasma 
potential  above  that  generated  with  a  purely 
absorbing  plate .  Also  presented  Is  a  parametric 
study  of  the  effects  of  secondary  electron  emission 
and  of  charged  particle  reflection  on  the  potential 
profile  and  on  particle  and  energy  transport  to  the 
wall . 


*  Work  supported  by  DOE  contract  DE-AT03-76ET53064 
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